Sleep Stage Detection
Systems App - (SSDApp)
Proyecto "Automatic Sleep Staging"
El presente notebook constituye el proyecto de detección y clasificación de las etapas de sueño por las que pasa un sujeto, de la asignatura Analítica de Datos en Salud, del grado en Ciencia de Datos de la Universtiat de València. Dicho proyecto, titulado "Sleep Stage Detection Systems App" y realizado por José Javier Gutiérrez Gil y Carmen Guarner Giner, se enmarca en el compromiso de explorar soluciones avanzadas en el campo de la salud, específicamente en el análisis automático de los estados del sueño.
El objetivo principal del proyecto es desarrollar un modelo de aprendizaje automático (Machine Learning, ML) que tenga la capacidad de clasificar los diferentes estados de sueño a través del análisis de un conjunto de señales provenientes del polisomnograma de un sujeto. Con el fin de cumplir dicho objetivo general, los objetivos específicos que determinamos son los siguientes.
Objetivo 1. Encontrar la mejor combinación de señales del polisomnograma dado; esto es, la combinación más adecuada para conseguir dicho objetivo. Notar que para todos los sujetos (subjects) de estudio se tendrá el mismo conjunto de señales.
Objetivo 2. Estudio en tiempo y frecuencia de las señales elegidas, con el fin de analizar la adecucación de aplicación de filtrados u otras operaciones sobre las mismas.
Objetivo 3. Extracción de características de las señales seleccionadas. Cálculo de un vector de características. Obtenión de características por trasnformadas (tiempo, frecuencia, tiempo-frecuencia) y aplicación de una normalización posterior de los resultados.
Objetivo 4. Entrenamiento del modelo ML de clasificación con los datos de los sujetos de estudio.
Objetivo 5. Evaluación de la calidad del modelo con el conjunto de sujetos de prueba.
Objetivo 6. Análisis de resultados, deliberación y concreción de conclusiones sobre los resultados obtenidos en las predicciones del modelo final.
El conjunto de datos con el que se trabaja para crear la aplicación y entrenar el modelo es el almacenado en el banco de datos ISRUC-Sleep. En concreto, se utiliza la serie Subgroup_3.
Los datos del conjunto de pacientes Subgroup_3 se presentan en carpetas comprimidas, cada uno compuesto por varios archivos. En primer lugar, encontramos un archivo .rec que es, en realidad, un archivo .edf, que debe ser renombrado para su correcta identificación. Además, cada carpeta comprimida incluye dos archivos .txt que contienen las anotaciones de los especialistas. Estos archivos son relevantes para la interpretación de las marcas asociadas a las distintas etapas del sueño. Asimismo, se proporcionan dos archivos .xlsx que contienen la misma información que los archivos .txt, pero también incorporan datos adicionales. Estos archivos .xlsx resultan especialmente útiles para realizar un análisis más profundo, permitiendo indagar en las razones detrás de posibles errores, descartar épocas de calidad dudosa, entre otros aspectos.
La serie Subgroup_3 consta de 10 sujetos (9 hombres y 1 mujer), con una sesión de adquisición de datos por sujeto. Se trata de sujetos saludables (grupo de control). Los años de los sujetos se sitúan entre 30 y 58 años, con el promedio ubicándose en 40 años y con una desviación estándar de 10 años.
El polisomnograma (PSG) asociado a estos datos incluye una variedad de señales provenientes de diferentes canales. Estos canales son:
Específicamente, cada grabación consta de señales de 19 canales. Todas las señales de EEG, EOG y EMG de mentón se muestrean a 200 Hz y se almacenan utilizando los formatos de datos estándar .edf. Además, todos los registros del conjunto de datos se encuentran segmentados en épocas de 30 segundos y puntuados visualmente por dos expertos en sueño diferentes en CHUC de acuerdo con las pautas de AASM, con las etapas:
$$Despierto - NREM (N1, N2 \text{ y } N3) - REM$$El registro para cada sujeto es el siguiente:
| Channel number | Type of the signal | Label | Frequency rate/Hz | Butterworth | Notch filter | Description |
|---|---|---|---|---|---|---|
| 1 | EOG | LOC-A2 | 200 | 0.3 Hz–35 Hz | 50 Hz | Left eyes MOV |
| 2 | ROC-A1 | Right eyes MOV | ||||
| 3 | F3-A2 | 50 Hz | Brain channels | |||
| 4 | C3-A2 | with the left | ||||
| 5 | EEG | O1-A2 | 200 | 0.3 Hz–35 Hz | 50 Hz | |
| 6 | F4-A1 | |||||
| 7 | C4-A1 | |||||
| 8 | O2-A1 | |||||
| 9 | Chin EMG | X1 | 200 | 10 Hz–70 Hz | 50 Hz | Chin EMG, placed |
| 10 | ECG (EKG) | X2 | 200 | Electrocardiographic | ||
| 11 | Leg-1 EMG | X3 | 200 | Left leg movement | ||
| 12 | Leg-2 EMG | X4 | 200 | 10 Hz–70 Hz | 50 Hz | Right leg movement |
| 13 | Snore | X5 | 200 | 10 Hz–70 Hz | 50 Hz | Snore (derived) |
| 14 | Flow-1 | X6 | 12.5 | |||
| 15 | Flow-2 | DC3 | 25 | Airflow (pressure based) | ||
| 16 | Abdominal | X7 | 25 | Abdominal efforts | ||
| 17 | Abdominal | X8 | 25 | Abdominal efforts | ||
| 18 | Pulse oximetry | SaO2 | 12.5 | Pulse oximetry (SaO2) | ||
| 19 | Body position | DC8 | 25 | Body position (BPOS) |
Las etapas del sueño se presentan a continuación, junto con las características asociadas a cada una de estas.
Las ondas que podemos encontrar en cada una de las etapas del sueño son:
| Ondas Estados de Sueño |
|---|
| Bandas | F (Hz) |
|---|---|
| low delta | 0.3–1 Hz |
| delta | 1–4 Hz |
| theta | 4–8 Hz |
| alpha | 8–12 Hz |
| sigma | 12–15 Hz |
| beta | 15–30 Hz |
| Fase del Sueño | Características | Señales del Polisomnograma | Ondas Asociadas |
|---|---|---|---|
| N1 | Ondas theta, breve | - EOG: LOC-A2 (Left eyes MOV), ROC-A1 (Right eyes MOV) | Theta (EEG) |
| - EEG: F3-A2, C3-A2 (Brain channels with the left) | |||
| N2 | Complejos K, husos, Theta, Delta | - EEG: O1-A2, F4-A1, C4-A1, O2-A1 | Theta, Delta (EEG) |
| N3 | Ondas delta, sueño profundo | - Chin EMG: X1 (Chin EMG, placed) | Delta (EEG) |
| REM | Movimientos oculares rápidos, actividad cerebral similar a la vigilia | - EEG: O1-A2, F4-A1, C4-A1, O2-A1 | Beta (EEG), Theta (EOG) |
| - EOG: LOC-A2, ROC-A1 | |||
| - EMG: X1 (Chin EMG, placed) | |||
| - ECG (EKG): X2 (Electrocardiographic) | |||
| - Leg-1 EMG: X3 (Left leg movement) | |||
| - Leg-2 EMG: X4 (Right leg movement) | |||
| - Snore: X5 (Snore, derived) | |||
| - Flow-1: X6, Flow-2: DC3 (Airflow, pressure based) | |||
| - Abdominal: X7, X8 (Abdominal efforts) | |||
| - Pulse oximetry: SaO2 (Pulse oximetry, SaO2) | |||
| - Body position: DC8 (Body position, BPOS) |
En la actual subsección, exponemos los patrones que podemos encontrar y distinguir una vez detectado las etapas de sueño en la señal del sujeto, así como distintas afecciones del tipo de sueño y el efecto de la medicación en las etapas del sueño.
Las afecciones de tipo de sueño se presentan seguidamente, diferenciando según la calidad, la duración y el patrón del sueño.
Categorías de Afecciones:
Insomnio:
Trastornos del Sueño del Ritmo Circadiano:
Parasomnias:
Trastornos Respiratorios del Sueño:
Narcolepsia y Trastornos Relacionados:
Apnea del sueño:
| Fase | Características | Ondas Cerebrales | Afecciones |
|---|---|---|---|
| N1 | Ondas theta, breve | Theta | Posible inicio de insomnio |
| N2 | Complejos K, husos del sueño | Theta, Delta | Insomnio de mantenimiento, trastornos del |
| sueño del ritmo circadiano | |||
| N3 | Ondas delta, sueño profundo | Delta | Ronquidos, apnea del sueño |
| REM | Movimientos oculares rápidos | Patrón similar a | Narcolepsia, trastornos del sueño del |
| Actividad cerebral similar a | vigilia, beta | comportamiento REM | |
| la vigilia |
Por último, queremos indicidir en los efectos de la medicación sobre las etapas del sueño. El uso crónico de benzodiazepinas induce un aumento de la etapa N2 del sueño, una disminución de la N3 (menor actividad delta y theta) y un aumento de los despertares. Del mismo modo, se puede observar la reducción del sueño REM o incluso una supresión completa en humanos después de la administración de antidepresivos tricíclicos. Además, los antidepresivos tenorarenérgicos y serotoninérgicos específicos mirtazapinas aumentan el tiempo total de sueño.
El presente notebook ha sido ejecutado desde Google Colab. Es por ello que incluimos algunos bloques de código para su configuración en dicha plataforma.
from google.colab import drive drive.mount('/content/drive')
import os os.chdir('/content/drive/MyDrive') os.chdir('/content/drive/MyDrive/base')
import sys sys.path.append('/content/drive/MyDrive')
# Instala la librería utilizando !pip
!pip install mne
!pip install yasa
!pip install dotmap
!pip install numpy
!pip install PyWavelets
!pip install statsmodels
!pip install scikit-learn==1.2.0
!pip install --upgrade scipy==1.8.1
#!pip install --upgrade pandas
Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: mne in /home/javi/.local/lib/python3.10/site-packages (1.6.1) Requirement already satisfied: scipy>=1.7.1 in /home/javi/.local/lib/python3.10/site-packages (from mne) (1.8.1) Requirement already satisfied: numpy>=1.21.2 in /home/javi/.local/lib/python3.10/site-packages (from mne) (1.24.4) Requirement already satisfied: decorator in /usr/lib/python3/dist-packages (from mne) (4.4.2) Requirement already satisfied: tqdm in /home/javi/.local/lib/python3.10/site-packages (from mne) (4.66.1) Requirement already satisfied: jinja2 in /usr/lib/python3/dist-packages (from mne) (3.0.3) Requirement already satisfied: lazy-loader>=0.3 in /home/javi/.local/lib/python3.10/site-packages (from mne) (0.3) Requirement already satisfied: matplotlib>=3.5.0 in /usr/lib/python3/dist-packages (from mne) (3.5.1) Requirement already satisfied: pooch>=1.5 in /home/javi/.local/lib/python3.10/site-packages (from mne) (1.8.0) Requirement already satisfied: packaging in /usr/lib/python3/dist-packages (from mne) (21.3) Requirement already satisfied: platformdirs>=2.5.0 in /home/javi/.local/lib/python3.10/site-packages (from pooch>=1.5->mne) (4.1.0) Requirement already satisfied: requests>=2.19.0 in /home/javi/.local/lib/python3.10/site-packages (from pooch>=1.5->mne) (2.31.0) Requirement already satisfied: certifi>=2017.4.17 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.19.0->pooch>=1.5->mne) (2023.11.17) Requirement already satisfied: urllib3<3,>=1.21.1 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.19.0->pooch>=1.5->mne) (2.1.0) Requirement already satisfied: charset-normalizer<4,>=2 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.19.0->pooch>=1.5->mne) (3.3.2) Requirement already satisfied: idna<4,>=2.5 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.19.0->pooch>=1.5->mne) (3.6) Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: yasa in /home/javi/.local/lib/python3.10/site-packages (0.6.4) Requirement already satisfied: pandas in /home/javi/.local/lib/python3.10/site-packages (from yasa) (2.1.4) Requirement already satisfied: mne>=1.3 in /home/javi/.local/lib/python3.10/site-packages (from yasa) (1.6.1) Requirement already satisfied: scikit-learn in /home/javi/.local/lib/python3.10/site-packages (from yasa) (1.2.0) Requirement already satisfied: numba>=0.57.1 in /home/javi/.local/lib/python3.10/site-packages (from yasa) (0.58.1) Requirement already satisfied: ipywidgets in /usr/lib/python3/dist-packages (from yasa) (6.0.0) Requirement already satisfied: lightgbm in /home/javi/.local/lib/python3.10/site-packages (from yasa) (4.2.0) Requirement already satisfied: scipy in /home/javi/.local/lib/python3.10/site-packages (from yasa) (1.8.1) Requirement already satisfied: sleepecg>=0.5.0 in /home/javi/.local/lib/python3.10/site-packages (from yasa) (0.5.6) Requirement already satisfied: seaborn in /home/javi/.local/lib/python3.10/site-packages (from yasa) (0.13.1) Requirement already satisfied: joblib in /home/javi/.local/lib/python3.10/site-packages (from yasa) (1.3.2) Requirement already satisfied: tensorpac>=0.6.5 in /home/javi/.local/lib/python3.10/site-packages (from yasa) (0.6.5) Requirement already satisfied: antropy in /home/javi/.local/lib/python3.10/site-packages (from yasa) (0.1.6) Requirement already satisfied: lspopt in /home/javi/.local/lib/python3.10/site-packages (from yasa) (1.3.0) Requirement already satisfied: numpy>=1.18.1 in /home/javi/.local/lib/python3.10/site-packages (from yasa) (1.24.4) Requirement already satisfied: matplotlib in /usr/lib/python3/dist-packages (from yasa) (3.5.1) Requirement already satisfied: pyriemann>=0.2.7 in /home/javi/.local/lib/python3.10/site-packages (from yasa) (0.5) Requirement already satisfied: tqdm in /home/javi/.local/lib/python3.10/site-packages (from mne>=1.3->yasa) (4.66.1) Requirement already satisfied: pooch>=1.5 in /home/javi/.local/lib/python3.10/site-packages (from mne>=1.3->yasa) (1.8.0) Requirement already satisfied: packaging in /usr/lib/python3/dist-packages (from mne>=1.3->yasa) (21.3) Requirement already satisfied: lazy-loader>=0.3 in /home/javi/.local/lib/python3.10/site-packages (from mne>=1.3->yasa) (0.3) Requirement already satisfied: jinja2 in /usr/lib/python3/dist-packages (from mne>=1.3->yasa) (3.0.3) Requirement already satisfied: decorator in /usr/lib/python3/dist-packages (from mne>=1.3->yasa) (4.4.2) Requirement already satisfied: llvmlite<0.42,>=0.41.0dev0 in /home/javi/.local/lib/python3.10/site-packages (from numba>=0.57.1->yasa) (0.41.1) Requirement already satisfied: PyYAML>=5.4.0 in /usr/lib/python3/dist-packages (from sleepecg>=0.5.0->yasa) (5.4.1) Requirement already satisfied: requests>=2.25.0 in /home/javi/.local/lib/python3.10/site-packages (from sleepecg>=0.5.0->yasa) (2.31.0) Requirement already satisfied: stochastic in /home/javi/.local/lib/python3.10/site-packages (from antropy->yasa) (0.7.0) Requirement already satisfied: pytz>=2020.1 in /usr/lib/python3/dist-packages (from pandas->yasa) (2022.1) Requirement already satisfied: tzdata>=2022.1 in /home/javi/.local/lib/python3.10/site-packages (from pandas->yasa) (2023.4) Requirement already satisfied: python-dateutil>=2.8.2 in /home/javi/.local/lib/python3.10/site-packages (from pandas->yasa) (2.8.2) Requirement already satisfied: threadpoolctl>=2.0.0 in /home/javi/.local/lib/python3.10/site-packages (from scikit-learn->yasa) (3.2.0) Requirement already satisfied: platformdirs>=2.5.0 in /home/javi/.local/lib/python3.10/site-packages (from pooch>=1.5->mne>=1.3->yasa) (4.1.0) Requirement already satisfied: six>=1.5 in /usr/lib/python3/dist-packages (from python-dateutil>=2.8.2->pandas->yasa) (1.16.0) Requirement already satisfied: charset-normalizer<4,>=2 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.25.0->sleepecg>=0.5.0->yasa) (3.3.2) Requirement already satisfied: urllib3<3,>=1.21.1 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.25.0->sleepecg>=0.5.0->yasa) (2.1.0) Requirement already satisfied: certifi>=2017.4.17 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.25.0->sleepecg>=0.5.0->yasa) (2023.11.17) Requirement already satisfied: idna<4,>=2.5 in /home/javi/.local/lib/python3.10/site-packages (from requests>=2.25.0->sleepecg>=0.5.0->yasa) (3.6) Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: dotmap in /home/javi/.local/lib/python3.10/site-packages (1.3.30) Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: numpy in /home/javi/.local/lib/python3.10/site-packages (1.24.4) Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: PyWavelets in /home/javi/.local/lib/python3.10/site-packages (1.5.0) Requirement already satisfied: numpy<2.0,>=1.22.4 in /home/javi/.local/lib/python3.10/site-packages (from PyWavelets) (1.24.4) Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: statsmodels in /home/javi/.local/lib/python3.10/site-packages (0.14.1) Requirement already satisfied: pandas!=2.1.0,>=1.0 in /home/javi/.local/lib/python3.10/site-packages (from statsmodels) (2.1.4) Requirement already satisfied: packaging>=21.3 in /usr/lib/python3/dist-packages (from statsmodels) (21.3) Requirement already satisfied: scipy!=1.9.2,>=1.4 in /home/javi/.local/lib/python3.10/site-packages (from statsmodels) (1.8.1) Requirement already satisfied: patsy>=0.5.4 in /home/javi/.local/lib/python3.10/site-packages (from statsmodels) (0.5.6) Requirement already satisfied: numpy<2,>=1.18 in /home/javi/.local/lib/python3.10/site-packages (from statsmodels) (1.24.4) Requirement already satisfied: tzdata>=2022.1 in /home/javi/.local/lib/python3.10/site-packages (from pandas!=2.1.0,>=1.0->statsmodels) (2023.4) Requirement already satisfied: python-dateutil>=2.8.2 in /home/javi/.local/lib/python3.10/site-packages (from pandas!=2.1.0,>=1.0->statsmodels) (2.8.2) Requirement already satisfied: pytz>=2020.1 in /usr/lib/python3/dist-packages (from pandas!=2.1.0,>=1.0->statsmodels) (2022.1) Requirement already satisfied: six in /usr/lib/python3/dist-packages (from patsy>=0.5.4->statsmodels) (1.16.0) Defaulting to user installation because normal site-packages is not writeable Requirement already satisfied: scikit-learn==1.2.0 in /home/javi/.local/lib/python3.10/site-packages (1.2.0) Requirement already satisfied: numpy>=1.17.3 in /home/javi/.local/lib/python3.10/site-packages (from scikit-learn==1.2.0) (1.24.4) Requirement already satisfied: threadpoolctl>=2.0.0 in /home/javi/.local/lib/python3.10/site-packages (from scikit-learn==1.2.0) (3.2.0) Requirement already satisfied: scipy>=1.3.2 in /home/javi/.local/lib/python3.10/site-packages (from scikit-learn==1.2.0) (1.8.1) Requirement already satisfied: joblib>=1.1.1 in /home/javi/.local/lib/python3.10/site-packages (from scikit-learn==1.2.0) (1.3.2)
import os
import sys
import warnings
import pandas as pd
import time
import numpy as np
import matplotlib.pyplot as plt
import mne
import yasa
from mne.datasets.brainstorm import bst_auditory
from mne.io import read_raw_ctf
from mne.preprocessing import annotate_movement, compute_average_dev_head_t
import numpy as np
from scipy.signal import resample
from scipy.signal import welch
Seguidamente, incluimos las clases definidas.
Config. Carga fichero de confguración ('config.yalm').from base.Config import Config
ExtractFeatures. Clase que contiene las funciones necesarias para extraer las características de las señal EEG y EOG indicadas en el artículo de referencia [1].from base.ExractFeatures import ExractFeatures
Polysom. Envuelve los datos del PSG leído de un sujeto determinado. Contiene funciones básicas de acceso y manipulación de las señales de cada uno de los canales del PSG. Se lee el fichero y se serializa el objeto mne para no mantenerlo en memoria RAM. Sólo se deserializa cuando se necesitan los datos. Únicamente se mantienen en memoria los datos de los canales eeg (eeg_raw_resam), los datos de los canales eog_eeg (eog_eeg_raw_resam) y la matriz de carateristicas normalizada (feature_norm_matrix).from base.Polysom import Polysom
Polysom_Data. Carga todos los ficheros de datos. Cada fichero contiene el polisomnograma de un sujeto determinado. Mantiene en memoria una lista de objetos Polysom, uno por cada sujeto (subject).try:
from base.Polysom_Data import Polysom_Data
except Exception as e:
print(f"Error importing module: {e}")
!ls /content/drive/MyDrive/base/config_psd.yalm
!ls /content/drive/MyDrive/base/data
!ls /content/drive/MyDrive/base/data/1
/content/drive/MyDrive/base/config_psd.yalm 1 10 2 3 4 5 6 7 8 9 ISRUC 1_1.txt 1_1.xlsx 1_2.txt 1_2.xlsx 1.edf
DEBUG = False
# Cargamos los polisomnogramas
pp = Polysom_Data()
iPath: ./data d: ./data/6 Polysom (d, subject, self.conf_data) file_path: ./data/6 subject: 6 path: ./data/6/6.edf Extracting EDF parameters from /home/javi/py1/data/6/6.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 5117999 = 0.000 ... 25589.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 8.8s
d: ./data/4 Polysom (d, subject, self.conf_data) file_path: ./data/4 subject: 4 path: ./data/4/4.edf Extracting EDF parameters from /home/javi/py1/data/4/4.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 4763999 = 0.000 ... 23819.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 8.4s
d: ./data/ISRUC invalid literal for int() with base 10: 'ISRUC' d: ./data/8 Polysom (d, subject, self.conf_data) file_path: ./data/8 subject: 8 path: ./data/8/8.edf Extracting EDF parameters from /home/javi/py1/data/8/8.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 5999999 = 0.000 ... 29999.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 10.2s
d: ./data/3 Polysom (d, subject, self.conf_data) file_path: ./data/3 subject: 3 path: ./data/3/3.edf Extracting EDF parameters from /home/javi/py1/data/3/3.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 4943999 = 0.000 ... 24719.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 8.5s
d: ./data/config_psd.yalm d: ./data/7 Polysom (d, subject, self.conf_data) file_path: ./data/7 subject: 7 path: ./data/7/7.edf Extracting EDF parameters from /home/javi/py1/data/7/7.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 4883999 = 0.000 ... 24419.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 8.5s
d: ./data/5 Polysom (d, subject, self.conf_data) file_path: ./data/5 subject: 5 path: ./data/5/5.edf Extracting EDF parameters from /home/javi/py1/data/5/5.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 5663999 = 0.000 ... 28319.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 9.8s
d: ./data/2 Polysom (d, subject, self.conf_data) file_path: ./data/2 subject: 2 path: ./data/2/2.edf Extracting EDF parameters from /home/javi/py1/data/2/2.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 5645999 = 0.000 ... 28229.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 9.6s
d: ./data/9 Polysom (d, subject, self.conf_data) file_path: ./data/9 subject: 9 path: ./data/9/9.edf Extracting EDF parameters from /home/javi/py1/data/9/9.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 5813999 = 0.000 ... 29069.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 9.9s
d: ./data/1 Polysom (d, subject, self.conf_data) file_path: ./data/1 subject: 1 path: ./data/1/1.edf Extracting EDF parameters from /home/javi/py1/data/1/1.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 5723999 = 0.000 ... 28619.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 9.6s
d: ./data/10 Polysom (d, subject, self.conf_data) file_path: ./data/10 subject: 10 path: ./data/10/10.edf Extracting EDF parameters from /home/javi/py1/data/10/10.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 4775999 = 0.000 ... 23879.995 secs... Filtering raw data in 1 contiguous segment Setting up low-pass filter at 33 Hz FIR filter parameters --------------------- Designing a one-pass, zero-phase, non-causal lowpass filter: - Windowed time-domain design (firwin) method - Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation - Upper passband edge: 33.33 Hz - Upper transition bandwidth: 8.33 Hz (-6 dB cutoff frequency: 37.50 Hz) - Filter length: 81 samples (0.405 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 8.5s
for sub in pp.lSubject:
print('Subject:', sub.subject)
sub.plot_signal_seg(seg = 20)
Subject: 6
Subject: 4
Subject: 8
Subject: 3
Subject: 7
Subject: 5
Subject: 2
Subject: 9
Subject: 1
Subject: 10
%matplotlib inline
for sub in pp.lSubject:
# sub.plot_signal_psd(average = True)
sub.plot_signal_psd()
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:495: UserWarning: Zero value in spectrum for channel DC8 spectrum = raw.compute_psd () /home/javi/py1/base/Polysom.py:497: UserWarning: Infinite value in PSD for channel DC8. These channels might be dead. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:495: UserWarning: Zero value in spectrum for channel DC8 spectrum = raw.compute_psd () /home/javi/py1/base/Polysom.py:497: UserWarning: Infinite value in PSD for channel DC8. These channels might be dead. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Effective window size : 10.240 (s)
/home/javi/py1/base/Polysom.py:497: RuntimeWarning: Channel locations not available. Disabling spatial colors. fig = spectrum.plot(average=average, picks="data", exclude="bads") /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
Indicaciones: Observamos que en el dominio del tiempo existen señales con artefactos que distorsinaran los valores. De hecho, esa parte de la señal no debería usarse para el estudio. Por otro lado, vemos en el dominio de la frecuencia que los archivos 3 y 8 tienen un PSD constante de casi -200. Esto se debe a que los canales DC8 Body position (BPOS) y DC3 son cero desde el inicio de la captación de la señal.
path = "data/%s/%s.edf" % (3, 3)
raw_8 = mne.io.read_raw_edf(path, preload = True)
spectrum = raw_8.compute_psd()
spectrum = raw_8.compute_psd(fmin = 0.5,
fmax = 30,
picks = "data",
n_fft = 2048)
spectrum.plot(average = True)
plt.show()
Extracting EDF parameters from /home/javi/py1/data/3/3.edf... EDF file detected Setting channel info structure... Creating raw.info structure... Reading 0 ... 4943999 = 0.000 ... 24719.995 secs... Effective window size : 10.240 (s)
/tmp/ipykernel_325790/298652374.py:5: UserWarning: Zero value in spectrum for channel DC8 spectrum = raw_8.compute_psd ()
Effective window size : 10.240 (s)
/tmp/ipykernel_325790/298652374.py:6: UserWarning: Zero value in spectrum for channel DC8 spectrum = raw_8.compute_psd (fmin=0.5, fmax=30, picks="data", n_fft=2048) /tmp/ipykernel_325790/298652374.py:7: UserWarning: Infinite value in PSD for channel DC8. These channels might be dead. spectrum.plot (average=True) /home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
mean_data = raw_8.get_data().mean(axis = 1)
std_data = raw_8.get_data().std(axis = 1)
import matplotlib.pyplot as plt
plt.hist(raw_8.get_data().flatten(),
bins = 100)
plt.xlabel('Valor')
plt.ylabel('Frecuencia')
plt.show()
mean_data
array([-1.34358095e-10, -2.14883293e-09, 3.26710659e-10, 6.06299674e-10,
-1.84286794e-10, 4.61806242e-10, -1.44365872e-09, -1.76396901e-09,
1.64442462e-09, -1.54694909e-08, 6.31256458e-10, 4.76613143e-09,
1.50878044e-08, -1.71373660e-08, 3.17630634e-05, -6.64014408e-08,
-1.92855641e-07, 9.54218784e+01, 0.00000000e+00])
spectrum = raw_8.copy().pick(['ROC-A1',
'LOC-A2',
'X1',
'X3',
'X4',
'X2']).compute_psd(fmin = 0.3,
fmax = 50)
spectrum.plot(average = True)
Effective window size : 10.240 (s)
/home/javi/.local/lib/python3.10/site-packages/mne/viz/utils.py:165: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. (fig or plt).show(**kwargs)
spectrum = raw_8.copy().pick(['O1-A2',
'F4-A1',
'C4-A1',
'O2-A1',
'F3-A2',
'C3-A2']).compute_psd(fmin = 0.3,
fmax = 50)
spectrum.plot()
Effective window size : 10.240 (s)
<ipython-input-19-53cf71e855aa>:2: RuntimeWarning: Channel locations not available. Disabling spatial colors. spectrum.plot()
# 'DC8'
spectrum = raw_8.copy().pick(['F3-A2',
'C3-A2',
'O1-A2',
'F4-A1',
'C4-A1',
'O2-A1',
'ROC-A1',
'LOC-A2',
'X1',
'X2',
'X3',
'X4',
'X7',
'X8',
'SaO2',
'DC3']).compute_psd(fmin = 0.3,
fmax = 50)
spectrum.plot(average = True)
Effective window size : 10.240 (s)
Al cargar los polisomnogramas de cada sujeto (subject) realizamos un remuestreo de la señal. Sin embargo, este remuestreo se realiza en dos fases para evitar perder información relevante antes de segmentar la señal en épocas de 30 segundos. Por tanto, en vez de utilizar raw.resample(), como procedemos es como sigue:
Aplicamos un filtro paso bajo de los datos en $\frac{1}{3}$ por debajo de la frecuencia de muestreo deseada. Tomamos 100 Hz.
Diezmamos los datos después de la creación de época y pasamos el parámetro decim calculado al cargar el fichero al constructor.
También sabemos que la frecuencia de muestreo es de 200Hz y la bajamos a 100Hz.
En primer lugar, según el teorema de muestreo de Nyquist-Shannon, para evitar el aliasing, la frecuencia de corte de un filtro paso bajo debe ser menor o igual a la mitad de la nueva frecuencia de muestreo. En este caso, la nueva frecuencia de muestreo es de 100 Hz, por lo que la frecuencia de corte no debe ser mayor de 50 Hz. Elegir una frecuencia de corte en 1/3 de la nueva frecuencia de muestreo (33.3 Hz) cumple con esta condición.
Así mismo, dado que la señal beta es de 30 Hz, un filtro paso bajo con una frecuencia de corte de 33.3 Hz aún permite que las componentes de la señal de interés pasen sin atenuación significativa. Esto significa que se conserva la información importante en los datos, mientras se atenua eficazmente las frecuencias más altas que pueden considerarse ruido o componentes no deseadas.
Además, reducir la frecuencia de muestreo a 100 Hz en lugar de mantenerla en 200 Hz disminuye la cantidad de datos que se debe procesar, lo que puede hacer que el análisis sea más eficiente sin sacrificar la información crítica.
import numpy as np
from sklearn.decomposition import FastICA, PCA
import matplotlib.pyplot as plt
n_subjects = len(pp.lSubject)
# psg_data' es una matriz de datos PSG (canales x muestras)
for i in range(n_subjects):
print(pp.lSubject[i].miscellaneous())
the (cropped) sample data object has 5118000 time samples and 19 channels. The last time sample is at 25589.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-03-28 00:10:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 4764000 time samples and 19 channels. The last time sample is at 23819.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-03-14 23:59:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 6000000 time samples and 19 channels. The last time sample is at 29999.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-02-28 00:00:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 4944000 time samples and 19 channels. The last time sample is at 24719.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-02-14 01:15:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 4884000 time samples and 19 channels. The last time sample is at 24419.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-03-14 23:16:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 5664000 time samples and 19 channels. The last time sample is at 28319.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-02-14 00:15:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 5646000 time samples and 19 channels. The last time sample is at 28229.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-04-25 00:21:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 5814000 time samples and 19 channels. The last time sample is at 29069.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-04-25 23:44:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 5724000 time samples and 19 channels. The last time sample is at 28619.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-03-28 23:19:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None the (cropped) sample data object has 4776000 time samples and 19 channels. The last time sample is at 23879.995 seconds. The first few channel names are LOC-A2, ROC-A1, F3-A2. bad channels: [] 200.0 Hz None <Info | 8 non-empty values bads: [] ch_names: LOC-A2, ROC-A1, F3-A2, C3-A2, O1-A2, F4-A1, C4-A1, O2-A1, X1, ... chs: 19 EEG custom_ref_applied: False highpass: 0.0 Hz lowpass: 33.3 Hz meas_date: 2009-03-07 00:32:00 UTC nchan: 19 projs: [] sfreq: 200.0 Hz subject_info: 1 item (dict) > None
La eliminación de artefactos del EEG comprende una serie de métodos para la reducción o cancelación de ruido, producido por señales subyacentes asociadas a procesos fisiológicos originadas por la misma fuente.
En el contexto de señales del cerebro obtenidas a partir de un polisomnograma (PSG), un artefacto generalmente se refiere a cualquier interferencia no deseada o distorsión en la señal que puede deberse a factores externos o internos. Estos artefactos pueden afectar la calidad de las señales de EEG (electroencefalograma), que son cruciales para el análisis del sueño en un PSG.
Algunas categorías comunes de artefactos en señales de EEG de un PSG son:
Movimientos Oculares (EOG): Movimientos rápidos de los ojos pueden introducir artefactos en las señales de EEG. Estos artefactos pueden deberse a movimientos oculares involuntarios durante el sueño.
Movimientos Musculares (EMG): Actividades musculares, especialmente en la región facial y de la mandíbula, pueden generar artefactos en las señales de EEG. Los cambios en la impedancia eléctrica debido a movimientos musculares pueden afectar la calidad de las señales.
Ruido Eléctrico y Ambiente: Interferencias eléctricas, como las generadas por equipos eléctricos, luces y otros dispositivos electrónicos, pueden introducir ruido en las señales de EEG.
Posición del Electrodo: La colocación incorrecta de electrodos o problemas con la fijación de electrodos pueden causar artefactos en las señales de EEG.
Fenómenos Fisiológicos Anómalos: Algunos fenómenos fisiológicos anómalos, como la presencia de sudor o movimientos involuntarios, pueden afectar las señales de EEG.
La detección y mitigación de artefactos en señales de EEG son aspectos críticos del procesamiento de datos de PSG. Métodos avanzados, como el uso de algoritmos de eliminación de artefactos, filtrado adaptativo y técnicas de rechazo de artefactos, son empleados para mejorar la calidad de las señales y obtener resultados más precisos en el análisis del sueño.
Es importante realizar un monitoreo cuidadoso durante la adquisición de datos y utilizar métodos de preprocesamiento adecuados para minimizar el impacto de los artefactos en las señales de EEG durante el análisis del sueño.
Cómo hemos indicado, la eliminación de artefactos del EEG comprende una serie de métodos para la reducción o cancelación de ruido, producido por señales subyacentes asociadas a procesos fisiológicos originadas por la misma fuente o externos.
Existen tres tipos de artefactos que no son cerebrales de origen fisiológico como son el parpadeo de ojos, la actividad cardiaca y la actividad muscular. Y hay un cuarto tipo correspondiente al equipamiento y los electrodos con origen técnico. Cabe resaltar que es posible que se presenten ambos tipos de artefactos durante una lectura de EEG. Más aún existen ciertos tipos de artefactos que se pueden encontrar de manera más frecuente bajo cierto tipo de actividad como el sueño.
Movimiento de los ojos y parpadeo: El movimiento de los ojos produce actividad eléctrica, cuando es lo suficiente fuerte puede ser registrada en el EEG. El EOG refleja la actividad realizada entre la córnea y la retina durante el movimiento ocular. La principal interferencia del EOG con la señal depende de la proximidad del electrodo y la dirección hacia la que se esté moviendo del ojo. El comportamiento repetitivo puede ser fácilmente identificado, el EOG puede ser confundido con actividad lenta del EEG como actividad theta o delta. El movimiento ocular no solo está presente durante fase de actividad consciente también puede estar presente durante el sueño.
Actividad muscular: Otro artefacto común corresponde a la actividad generada por la contracción muscular, medida en la superficie del cuerpo con Electromiografía EMG, este tipo de artefacto es encontrado en pacientes despiertos al realizar comportamientos como tragar saliva, hacer muecas, fruncir el ceño, masticar, hablar, chupar, e hipo (Sörnmo & Laguna, 2005). La manera general de la forma del EMG depende del grado de contracción del musculo. Una contracción tenue produce un tren de picos de baja amplitud, mientras que un aumento en la fuerza de la contracción disminuye la distancia entre picos de tal manera que el EMG exhibe el comportamiento de una señal continuamente variante.
Actividad cardiaca: La actividad cardiaca del corazón se ve reflejada por el electrocardiograma ECG, no obstante, la actividad cardiaca presente en los registros de EEG suele ser particularmente baja. (1 - 2, 20-100 𝜇 V respectivamente). Puede obstaculizar la actividad EEG en ciertos tipos de fisionómicas como sujetos robustos con cuellos cortos y gruesos, la forma de la onda de patrón regular ayuda a caracterizar los latidos del corazón en condiciones normales, siendo muy útil en la cancelación de artefactos. Aunque en algunos casos las formas de pico pueden confundirse durante ataques epilépticos, y también en la presencia de arritmias cardiacas exhibiendo una amplia variación en el intervalo entre latidos.
Electrodos y equipamiento: El artefacto “electrodo-pop” de su traducción al español, se encuentra basado en la manifestación del contacto de la corriente directa DC con la interacción del electrodo sobre la piel. Se puede presentar en cualquier tipo de señal bio eléctrica, y es manifestado por un cambio abrupto sobre la línea base seguido por una reincorporación gradual a la misma. En ocasiones puede ser confundido con ondas en pico. Otro tipo de artefacto que ser posible es mediante el cable que conecta el dispositivo, ya que si el cable no cuenta con la adecuada cobertura lo hace susceptible a ruido electromagnético causado por corrientes de dispositivos electrónicos cercanos, como resultado una interferencia de potencia de 50 a 60 Hz es concentrada en los electrodos contaminando la señal EEG. Finalmente, los artefactos relacionados con el equipamiento como ruido interno de los amplificadores y recorte de amplitud causado por los conversores análogo – digital con rango dinámico muy estrecho.
import numpy as np
import matplotlib.pyplot as plt
# Definimos las señales independientes (S)
np.random.seed(42) # Fijamos la semilla para reproducibilidad
S = np.random.rand(4, 1000) # 4 señales independientes de longitud 1000
# Definimos la matriz de mezcla (A)
A = np.random.rand(4, 4)
# Generamos las señales observadas (X)
X = np.dot(A, S)
# Aplicamos ICA para descomponer X en A y S
from sklearn.decomposition import FastICA
ica = FastICA(n_components=4)
S_ = ica.fit_transform(X.T) # Transponemos X ya que sklearn ICA trabaja con datos transpuestos
# Recuperamos las señales independientes estimadas
X_ = np.dot(S_, ica.mixing_.T)
# Graficamos las señales originales, mezcladas y recuperadas
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.title("Señales Originales (Independientes)")
plt.plot(S.T)
plt.subplot(3, 1, 2)
plt.title("Señales Mezcladas (Observadas)")
plt.plot(X.T)
plt.subplot(3, 1, 3)
plt.title("Señales Recuperadas (Estimadas)")
plt.plot(X_.T)
plt.tight_layout()
plt.show()
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn(
Este algoritmo realmente funciona con las señales EEG y observamos su correlación con el movimiento de parpadeos EOG. Ahora se le pasan todos los canales.
import numpy as np
from sklearn.decomposition import FastICA, PCA
import matplotlib.pyplot as plt
# psg_data' es una matriz de datos PSG (canales x muestras)
psg_data = pp.lSubject[8].getData()
# print(type(psg_data))
# print(psg_data.info)
# print(psg_data.get_data())
# psg_data = psg_data.get_data()
# Reducción de dimensionalidad con PCA
# Número de componentes principales a mantener
n_components_pca = 0.99
pca = PCA(n_components = n_components_pca)
psg_data_reduced = pca.fit_transform(psg_data.T)
# Inicializa el modelo ICA
# Número de componentes independientes a extraer
n_components_ica = 6
ica = FastICA(n_components = n_components_ica,
random_state = 0)
# Ajustamos el modelo a los datos PSG reducidos
ica.fit(psg_data_reduced)
# Obtiene las componentes independientes
# psg_data_ica contiene las señales EEG transformadas, donde cada columna representa una componente independiente.
psg_data_ica = ica.transform(psg_data_reduced)
# Inverse transform para obtener las señales separadas
# para transformar las componentes independientes nuevamente a las señales originales
separated_signals = ica.inverse_transform(psg_data_ica)
# Calcula la desviación estándar de cada señal separada
std_devs = np.std(separated_signals, axis = 1)
# Visualizamos las señales independientes
plt.figure(figsize = (12, 4))
plt.plot(psg_data_ica)
plt.title('Señales Independientes de ICA')
plt.xlabel('Muestras')
plt.ylabel('Amplitud')
plt.show()
# Visualizamos las señales independientes
plt.figure(figsize = (12, 4))
plt.plot(separated_signals)
plt.title('Señales separated_signals de ICA')
plt.xlabel('Muestras')
plt.ylabel('Amplitud')
plt.show()
# Histograma de amplitudes de las señales independientes
plt.figure(figsize = (8, 6))
plt.hist(psg_data_ica.flatten(),
bins = 50,
color = 'c',
edgecolor = 'k',
alpha = 0.7)
plt.title('Histograma de Amplitudes de Señales Independientes')
plt.xlabel('Amplitud')
plt.ylabel('Frecuencia')
plt.show()
# Establecemos un umbral para decidir qué canales mostrar
# Calcular el umbral como una fracción de la desviación estándar
threshold_factor = 3 # Ajustamos el factor según tus necesidades
threshold = threshold_factor * np.std(psg_data_ica)
# Trama las señales separadas de los canales significativos
num_channels = psg_data.shape[0]
for channel in range(num_channels):
print(channel)
print(std_devs[channel])
if np.abs(std_devs[channel]) > threshold and separated_signals.T.shape [0] > channel:
plt.figure()
plt.plot(separated_signals.T[channel], label = f'Channel {channel + 1}')
plt.title(f'Separated Signal for Channel {channel + 1}')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()
plt.close()
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:589: UserWarning: n_components is too large: it will be set to 2 warnings.warn(
0 9.194477345025852 1 9.555170634236081 2 9.873985172948117 3 10.12566815896821 4 10.306320513854054 5 10.424797566849314 6 10.490120866047207 7 10.505457189751933 8 10.471068335084627 9 10.390171387984253 10 10.270812486132424 11 10.12278657328736 12 9.95398232338113 13 9.770046023977056 14 9.576777256959907 15 9.381820315345227 16 9.193713680853588 17 9.019670366107551 18 8.864585102553242
import numpy as np
import matplotlib.pyplot as plt
import mne
from sklearn.decomposition import PCA, FastICA
# Función para detectar artefactos usando PCA e ICA
def detect_artifacts (raw):
# Obtenemos datos brutos
data = raw.get_data()
# Aplicamos PCA
n_components_pca = 0.99 # Conservamos el 99% de la varianza
pca = PCA(n_components = n_components_pca)
data_pca = pca.fit_transform(data.T)
# Aplicamos ICA
n_components_ica = min(data.shape[0], 15) # Número máximo de componentes ICA
ica = FastICA(n_components=n_components_ica, random_state=97)
data_ica = ica.fit_transform(data_pca.T)
# Detectamos artefactos usando umbral en las señales independientes
threshold = 5 # Ajustamos el umbral según tus necesidades
artifact_indices = np.where(np.abs(data_ica) > threshold)
return artifact_indices
# Función para realizar anotaciones en base a los artefactos detectados
def annotate_artifacts(raw, artifact_indices):
# Convertimos tiempos de anotación a índices de muestras
sample_indices = np.asarray(raw.times[artifact_indices] * raw.info['sfreq'], dtype=int)
# Creamos objeto de anotaciones
annotations = mne.Annotations(onset = raw.times[artifact_indices],
duration = np.zeros_like(artifact_indices),
description = 'Artifact')
# Añadimos anotaciones al objeto raw
raw.set_annotations(annotations)
# Iteramos sobre todos los polisomnogramas en la lista
for sub in pp.lSubject:
# Obtenemos el objeto raw
raw = sub.getEEG()
# Detectamos artefactos
artifact_indices = detect_artifacts(raw)
print(artifact_indices)
# Verificamos si se detectaron artefactos
if len(artifact_indices[0]) == 0:
print(f"No se detectaron artefactos en {sub.subject}")
continue
# Realizamos anotaciones
annotate_artifacts(raw, artifact_indices [0])
# Visualizamos polisomnograma con anotaciones
fig, ax = plt.subplots()
raw.plot(start = 0,
duration = 60,
n_channels = 5,
events = raw.annotations,
show = False)
ax.scatter(raw.times[artifact_indices],
raw._data[:, artifact_indices],
color = 'red',
marker = 'o',
label = 'Artifact')
ax.legend()
plt.title(f'Polisomnograma con Artefactos Detectados - {sub.subject}')
plt.show()
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:589: UserWarning: n_components is too large: it will be set to 5 warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:123: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 6
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 4
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 8
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:589: UserWarning: n_components is too large: it will be set to 5 warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 3
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:589: UserWarning: n_components is too large: it will be set to 5 warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 7
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:123: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 5
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:589: UserWarning: n_components is too large: it will be set to 5 warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:123: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 2
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:589: UserWarning: n_components is too large: it will be set to 5 warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:123: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 9
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 1
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:123: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations. warnings.warn(
(array([], dtype=int64), array([], dtype=int64)) No se detectaron artefactos en 10
import numpy as np
from sklearn.decomposition import FastICA, PCA
import matplotlib.pyplot as plt
# psg_data' es una matriz de datos PSG (canales x muestras)
psg_data = pp.lSubject [8].getData ()
# Reducción de dimensionalidad con PCA
n_components_pca = 19 # número de componentes principales que deseas mantener
pca = PCA (n_components = n_components_pca)
psg_data_reduced = pca.fit_transform (psg_data.T)
# Inicializa el modelo ICA
n_components_ica = 3 # Número de componentes independientes que deseas extraer
ica = FastICA (n_components = n_components_ica, random_state = 0)
# Ajusta el modelo a los datos PSG reducidos
ica.fit (psg_data_reduced)
# Obtiene las componentes independientes
# psg_data_ica contiene las señales EEG transformadas, donde cada columna representa una componente independiente.
psg_data_ica = ica.transform (psg_data_reduced)
# Inverse transform para obtener las señales separadas
# para transformar las componentes independientes nuevamente a las señales originales
separated_signals = ica.inverse_transform (psg_data_ica)
# Calcula la desviación estándar de cada señal separada
std_devs = np.std (separated_signals, axis = 1)
# Establece un umbral para decidir qué canales mostrar
threshold = 0.1
# Trama las señales separadas de los canales significativos
num_channels = psg_data.shape[0]
for channel in range (num_channels):
if std_devs [channel] > threshold:
plt.figure ()
plt.plot (separated_signals.T[channel], label = f'Channel {channel + 1}')
plt.title (f'Separated Signal for Channel {channel + 1}')
plt.xlabel ('Sample')
plt.ylabel ('Amplitude')
plt.legend ()
plt.grid (True)
plt.show ()
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn(
Probamos el método detección de parpadeos en señales EEG.
psg_data_ica = None
ica = None
psg_data_reduced = None
del psg_data_ica, pca, ica, psg_data_reduced, separated_signals
import gc
gc.collect()
16442
import numpy as np
from sklearn.decomposition import FastICA, PCA
import matplotlib.pyplot as plt
threshold_parpadeo = 6 # Umbral para detección de parpadeo
# Inicializamos el modelo ICA
# Número de componentes independientes a extraer
n_components_ica = 6
# Detección de parpadeo en cada componente independiente
parpadeo_detected = np.zeros ((n_components_ica,), dtype=bool)
# psg_data' es una matriz de datos PSG (canales x muestras)
psg_data = pp.lSubject[8].getEEG () # pp.lSubject [8].getData ()
print(psg_data)
import numpy as np
from sklearn.decomposition import FastICA, PCA
import matplotlib.pyplot as plt
threshold_parpadeo = 6 # Umbral para detección de parpadeo
# Inicializamos el modelo ICA
n_components_ica = 6 # Número de componentes independientes A extraer
# Detección de parpadeo en cada componente independiente
parpadeo_detected = np.zeros((n_components_ica, ),
dtype=bool)
# psg_data' es una matriz de datos PSG (canales x muestras)
psg_data = pp.lSubject[8].getEEG().get_data() # pp.lSubject[8].getData() //.getEEG ().get_data ()
# Reducción de dimensionalidad con PCA
n_components_pca = 4 # Número de componentes principales a mantener
pca = PCA(n_components = n_components_pca)
psg_data_reduced = pca.fit_transform(psg_data.T)
# Obtenemos señales independientes (modelo ICA)
ica = FastICA(n_components = n_components_ica,
random_state = 0)
# Ajustamos el modelo a los datos PSG reducidos
ica.fit(psg_data_reduced)
# Obtenemos las componentes independientes
# psg_data_ica contiene las señales EEG transformadas, donde cada columna representa una componente independiente.
psg_data_ica = ica.transform(psg_data_reduced)
for componente in range(n_components_ica):
if (psg_data_ica.shape[1] > componente):
signal = psg_data_ica[:, componente]
diff_signal = np.diff(signal)
parpadeo_count = 0
for i in range (1, len (diff_signal) - 1):
if diff_signal [i - 1] > 0 and diff_signal [i] < 0 and diff_signal [i + 1] > 0:
parpadeo_count += 1
if parpadeo_count >= threshold_parpadeo:
parpadeo_detected [componente] = True
# 'parpadeo_detected' contendrá información binaria donde 1 indica la presencia de parpadeo en cada canal (cada componente).
num_components = psg_data_ica.shape[1]
for componente in range(num_components):
plt.figure()
plt.plot(psg_data_ica[:, componente],
label = f'Componente {componente + 1}')
if parpadeo_detected[componente]:
# Marcamos puntos rojos donde se detecta el parpadeo
parpadeo_indices = np.where(diff_signal > 0)[0] # Indices donde ocurre el patrón de parpadeo
plt.scatter(parpadeo_indices, psg_data_ica[parpadeo_indices,
componente],
color = 'red',
label = 'Parpadeo')
plt.title(f'Componente {componente + 1} con Parpadeo (Detectado)')
else:
plt.title(f'Componente {componente + 1}')
plt.xlabel('Muestra')
plt.ylabel('Amplitud')
plt.legend()
plt.grid(True)
plt.show()
num_components = parpadeo_detected.shape [0]
for componente in range (num_components):
if parpadeo_detected [componente]:
print(f'Parpadeo detectado en la Componente {componente + 1}')
else:
print(f'No se detectó parpadeo en la Componente {componente + 1}')
# Creamos anotaciones con MNE
sfreq = 200
anotaciones = mne.Annotations(onset = parpadeo_indices / sfreq, # Convertir índices a segundos
duration = 1.0, # Duración de cada anotación en segundos
description = 'parpadeo')
# Creamos objeto Raw de MNE
info = mne.create_info(ch_names = [f'Ch{i+1}' for i in range(4)],
sfreq = sfreq,
ch_types = 'eeg')
nn_ch = psg_data_ica.shape[1]
raw = mne.io.RawArray(psg_data_ica.T[:, :nn_ch], info)
# Asignamos anotaciones al objeto Raw
raw.set_annotations(anotaciones)
# Visualizamos PSG con señales biológicas y anotaciones
# raw.plot(events=anotaciones, scalings='auto', start=0, duration=10, show_options=True)
# Ajustamos anotaciones para que estén dentro del rango de los datos
raw.set_annotations(anotaciones.delete(anotaciones.duration > raw.times[-1]))
# Visualizamos PSG con señales biológicas y anotaciones
events, _ = mne.events_from_annotations(raw)
raw.plot(events = events, show_options = True)
plt.show()
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn( /home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:589: UserWarning: n_components is too large: it will be set to 4 warnings.warn( /usr/lib/python3/dist-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc="best" can be slow with large amounts of data. fig.canvas.print_figure(bytes_io, **kw)
Parpadeo detectado en la Componente 1
Parpadeo detectado en la Componente 2
Parpadeo detectado en la Componente 3
Parpadeo detectado en la Componente 4
No se detectó parpadeo en la Componente 5
No se detectó parpadeo en la Componente 6
Creating RawArray with float64 data, n_channels=4, n_times=4
Range : 0 ... 3 = 0.000 ... 0.015 secs
Ready.
/tmp/ipykernel_325790/1178868782.py:87: RuntimeWarning: Omitted 2893373 annotation(s) that were outside data range. raw.set_annotations(anotaciones)
Using matplotlib as 2D backend.
/tmp/ipykernel_325790/1178868782.py:87: RuntimeWarning: Limited 4 annotation(s) that were expanding outside the data range. raw.set_annotations(anotaciones)
psg_data_ica = None
ica = None
psg_data = None
psg_data_reduced = None
del psg_data, psg_data_ica, ica, psg_data_reduced
import gc
gc.collect()
10645
import numpy as np
from sklearn.decomposition import FastICA, PCA
import matplotlib.pyplot as plt
# 'psg_data' es una matriz de datos PSG (canales x muestras)
psg_data = pp.lSubject[8].getData ()
# Reducción de dimensionalidad con PCA
n_components_pca = 19 # Número de componentes principales a mantener. Será el número de canales que tengo de origen
pca = PCA(n_components = n_components_pca)
psg_data_reduced = pca.fit_transform (psg_data.T)
# Inicializamos el modelo ICA
n_components_ica = 10 # Número de componentes independientes que deseamos extraer
ica = FastICA(n_components = n_components_ica, random_state = 0)
# Ajustamos el modelo a los datos PSG reducidos
ica.fit(psg_data_reduced)
# Obtienemos las componentes independientes
independent_components = ica.transform(psg_data_reduced)
# Inverse transform para obtener las señales separadas
separated_signals = ica.inverse_transform(independent_components)
# Detección de artefactos (ejemplo simple: umbral de amplitud)
threshold = 20 # Define un umbral de amplitud para detectar artefactos en las señales separadas
# Creamos una lista para almacenar los índices de artefactos por canal
artifacts_indices = []
for channel in range(psg_data.shape [0]):
signal = separated_signals.T[channel]
# mos los índices donde la amplitud excede el umbral
artifact_indices = np.where(np.abs (signal) > threshold)[0]
artifacts_indices.append(artifact_indices)
# Trama de las señales separadas con artefactos resaltados
for channel, artifact_indices in enumerate(artifacts_indices):
plt.figure()
plt.plot(separated_signals.T[channel], label = f'Channel {channel + 1}')
plt.scatter(artifact_indices, separated_signals.T[channel][artifact_indices], c = 'r', label = 'Artifacts')
plt.titl(f'Separated Signal for Channel {channel + 1} with Artifacts')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()
/home/javi/.local/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:542: FutureWarning: Starting in v1.3, whiten='unit-variance' will be used by default. warnings.warn(
/usr/lib/python3/dist-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc="best" can be slow with large amounts of data. fig.canvas.print_figure(bytes_io, **kw)
yasahypno = pp.lSubject[2].get_raw_eData()
raw_ = pp.lSubject[2].getData()
data = raw_ #sujto_0.getData()
sf = 200
art, zscores = yasa.art_detect(data, sf, window = 30, method = 'covar', threshold = 3)
art.shape, zscores.shape
WARNING:yasa:Flat channel(s) were found and removed in data.
((824,), (824,))
print(f'{art.sum()} / {art.size} epochs rejected.')
0 / 941 epochs rejected.
plt.plot(art);
plt.yticks([0, 1], labels = ['Good (0)', 'Art (1)']);
import seaborn as sns
sns.histplot(zscores)
plt.title('Histogram of z-scores')
plt.xlabel('Z-scores')
plt.ylabel('Density')
plt.axvline(3, color = 'r', label = 'Threshold')
plt.axvline(-3, color = 'r')
plt.legend(frameon = False);
plt.show()
from scipy.special import erf
threshold = 3
perc_expected_rejected = (1 - erf(threshold / np.sqrt(2))) * 100
print(f'{perc_expected_rejected:.2f}% of all epochs are expected to be rejected.')
0.27% of all epochs are expected to be rejected.
# Actual
(art.sum() / art.size) * 100
print(f'{(art.sum() / art.size) * 100:.2f}% of all epochs were actually rejected.')
0.00% of all epochs were actually rejected.
import numpy as np
import yasa
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.special import erf
# Supongamos que pp.lSubject es una lista de objetos Subject
def detect_and_plot_artifacts(data, sf, window = 30, method = 'covar', threshold = 3):
# Detectar artefactos
art, zscores = yasa.art_detect(data,
sf,
window = window,
method = method,
threshold = threshold)
# Imprimimos información sobre los artefactos
print(f'\t\t{art.sum()} / {art.size} epochs rejected.')
# Graficamos la presencia de artefactos
plt.plot(art)
plt.yticks([0, 1], labels=['Good (0)', 'Art (1)'])
plt.show()
# Histograma de z-scores
sns.histplot(zscores)
plt.title('Histogram of z-scores')
plt.xlabel('Z-scores')
plt.ylabel('Density')
plt.axvline(threshold, color = 'r', label = 'Threshold')
plt.axvline(-threshold, color = 'r')
plt.legend(frameon = False)
plt.show()
# Calculamos y mostramos el porcentaje esperado de rechazo
perc_expected_rejected = (1 - erf(threshold / np.sqrt(2))) * 100
print(f'\t\t {perc_expected_rejected:.2f}% of all epochs are expected to be rejected.')
# Calculamos y mostramos el porcentaje real de rechazo
perc_actual_rejected = (art.sum() / art.size) * 100
print(f'\t\t {perc_actual_rejected:.2f}% of all epochs were actually rejected.')
return art
# Iteramos sobre todos los sujetos y obtenemos los índices de artefactos para cada uno
for subject in pp.lSubject:
raw_data = subject.getEEG()
artifact_indices = detect_and_plot_artifacts(raw_data, sf = 200)
if len(artifact_indices) > 0 and raw_data is not None:
if raw_data.annotations is not None:
# Eliminamos anotaciones anteriores (si las hay)
for idx in range(len(raw_data.annotations)):
raw_data.annotations.delete(idx)
raw_data.annotations.append(artifact_indices / 200, np.zeros_like(artifact_indices), 'artifact')
0 / 954 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.00% of all epochs were actually rejected. 3 / 941 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.32% of all epochs were actually rejected. 0 / 824 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.00% of all epochs were actually rejected. 0 / 794 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.00% of all epochs were actually rejected. 0 / 944 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.00% of all epochs were actually rejected. 2 / 853 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.23% of all epochs were actually rejected. 0 / 814 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.00% of all epochs were actually rejected. 14 / 1000 epochs rejected.
0.27% of all epochs are expected to be rejected. 1.40% of all epochs were actually rejected. 0 / 969 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.00% of all epochs were actually rejected. 0 / 796 epochs rejected.
0.27% of all epochs are expected to be rejected. 0.00% of all epochs were actually rejected.
Se filtra para eliminar el ruido no deseado, con el objetivo de mejorar la calidad de la señal PSG y aumentar la SNR (relación señal-ruido). También poder remuestrear:
1) Un filtro para eliminar el ruido eléctrico de 50Hz. NR: Los datos cargados ya tienen eliminada esta frecuencia.
2) Un filtro Butterworth de paso de banda con un corte inferior de frecuencia de 0.3 Hz y un corte superior de frecuencia de 45 Hz para canales EEG y para el resto de canales las frecuencias necesarias. Al igual que el paso anterior, este filtrado ya está hecho en los datos cargados, para cada sujeto.
3) Notar que al cargar el PSG de cada sujeto ya se ha realizado un filtrado pasa baja a 33 Hz. Ya que deseamos remuestrear posteriormente la señal a 100 Hz, siendo 33 Hz (1/3) * 100 = 33 Hz de la frecuencia de muestreo deseada.
Según el artículo de referencia [1], a notch filter at 50 Hz (los datos ya vienen con el fitro realizado) a band-pass Butterworth filter with lower cut-off of 0.5 Hz and higher cut-off of 45 Hz. (Idem).
En esta subsección, procedemos a la segmentación en ventanas de 30 segundos y diezmado de la señal pasando el valor dicem calculado al cargar el PSG.
Diezmamos la señal por ventana. La función doSegmentResample realmente primero segmenta y luego hace el remuestreo indicado en la creacion del objeto.
Por otro lado, una vez hecha la segmentación y el remuestreo, tal como nos indica el paper, eliminar las 30 últimas componentes (muestras) de cada señal por cada canal.
for sub in pp.lSubject:
# print ('Subject:',sub.subject)
eeg_raw = sub.getEEG()
eeg_names = eeg_raw.info['ch_names']
# print ('\t',eeg_names)
time, eeg_resam = sub.doSegmentResample (eeg_raw)
# print ('\t',eeg_resam.shape)
eeg_resam = eeg_resam[:-30, :, :] # Eliminamos las 30 últimas porque lo indica el artículo de referencia
sub.set_eeg_raw_resam(eeg_resam, eeg_names)
# print('\t',eeg_resam.shape)
['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1']
# Los canales EEG y EOG:
for sub in pp.lSubject:
# print('Subject:', sub.subject)
eeg_eog_raw = sub.getChannels(eeg = True, eog = True)
eeg_eog_names = eeg_eog_raw.info['ch_names']
# print('\t',eeg_eog_names)
time, eeg_eog_resam = sub.doSegmentResample(eeg_eog_raw)
# print('\t',eeg_eog_resam.shape)
eeg_eog_resam = eeg_eog_resam[:-30,:,:] # Eliminamos las 30 últimas porque lo indica el artículo de referencia
# print('\t',eeg_eog_resam.shape)
sub.set_eeg_eog_raw_resam(eeg_eog_resam, eeg_eog_names)
# print('\t',eeg_resam.shape)
['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2'] ['O1-A2', 'F4-A1', 'C4-A1', 'O2-A1', 'ROC-A1', 'LOC-A2', 'F3-A2', 'C3-A2']
La Clase ExtractFeatures tiene todas las funciones. Podemos extraer:
La relación de dichas características con las fases del sueño son:
EEG and EOG channels - Entropía (Sallis, Renyi, Shannon):
Harmonic Parameters y Hjorth Parameters:
Relative Spectral Power:
Slow Wave Index (SWI):
Autoregressive Coefficients (order 3):
Percentiles (25, 50, 75), Skewness, Kurtosis:
EOG channels - Peak to Peak Amplitude:
EEG y Each sub-band - MODWT-based Features: Energy, Percentage of Energy, Mean and Standard deviation of each sub-band:
extract = ExractFeatures(pp.conf_data)
for sub in pp.lSubject:
# print('Subject:', sub.subject)
eeg_resam, eeg_names = sub.get_eeg_raw_resam()
eeg_eog_resam, eeg_eog_names = sub.get_eeg_eog_raw_resam()
sub.set_paper_feature_norm_matrix(extract.extract_paper_features(eeg_resam,
eeg_names,
eeg_eog_resam,
eeg_eog_names,
level = 3))
pp.lSubject[0].get_paper_feature_norm_matrix()
Subject: 6 Subject: 4 Subject: 8 Subject: 3 Subject: 7 Subject: 5 Subject: 2 Subject: 9 Subject: 1 Subject: 10
| Energy_s_nivel-1_Canal-O1-A2 | Energy_d_nivel-1_Canal-O1-A2 | Percentage_energy_s_nivel-1_Canal-O1-A2 | Percentage_energy_d_nivel-1_Canal-O1-A2 | Mean_s_nivel-1_Canal-O1-A2 | Mean_d_nivel-1_Canal-O1-A2 | Std_deviation_s_nivel-1_Canal-O1-A2 | Std_deviation_d_nivel-1_Canal-O1-A2 | Energy_s_nivel-1_Canal-F4-A1 | Energy_d_nivel-1_Canal-F4-A1 | ... | Percentil50_F3-A2 | Percentil50_C3-A2 | Percentil75_O1-A2 | Percentil75_F4-A1 | Percentil75_C4-A1 | Percentil75_O2-A1 | Percentil75_ROC-A1 | Percentil75_LOC-A2 | Percentil75_F3-A2 | Percentil75_C3-A2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 7.877584e-13 | 3.827433e-12 | 1.793814 | 64.805300 | 1.730601e-07 | 8.025774e-07 | 4.086438e-07 | 5.592207e-07 | 9.984503e-13 | 4.443003e-12 | ... | 8.340305e-08 | 1.397417e-07 | 1.894502e-06 | 1.976343e-06 | 1.736273e-06 | 1.762761e-06 | 2.058265e-06 | 2.604764e-06 | 2.199003e-06 | 1.866567e-06 |
| 1 | 2.884275e-12 | 4.378524e-12 | 51.534521 | 55.712371 | 7.477975e-07 | -8.668120e-07 | 4.023278e-07 | 5.858907e-07 | 2.768917e-12 | 1.397024e-11 | ... | 3.287122e-08 | 9.796863e-09 | 1.388383e-06 | 8.380086e-07 | 8.486788e-07 | 1.091260e-06 | 7.302970e-07 | 6.950321e-07 | 8.387663e-07 | 9.265201e-07 |
| 2 | 3.763747e-12 | 1.942999e-12 | 8.909356 | 88.111377 | 6.647579e-07 | -6.323346e-08 | 7.064233e-07 | 6.940831e-07 | 3.972120e-12 | 3.437004e-12 | ... | -8.448397e-09 | -2.335692e-08 | 9.158655e-07 | 9.438726e-07 | 8.577809e-07 | 7.830601e-07 | 8.799506e-07 | 6.437996e-07 | 8.252971e-07 | 7.934583e-07 |
| 3 | 9.245739e-14 | 6.595357e-12 | 0.357688 | 76.459332 | 2.757269e-08 | -1.040333e-06 | 1.495129e-07 | 7.526923e-07 | 2.773305e-13 | 3.842782e-12 | ... | -1.883530e-08 | -5.736141e-09 | 6.934272e-07 | 7.717873e-07 | 7.070651e-07 | 6.423942e-07 | 7.773066e-07 | 6.506189e-07 | 8.466400e-07 | 7.569879e-07 |
| 4 | 8.251030e-13 | 3.228449e-11 | 5.452811 | 96.955454 | 4.417207e-07 | -2.804984e-06 | 1.056342e-07 | 4.507639e-07 | 3.283667e-13 | 3.152991e-11 | ... | 3.462841e-08 | 2.123687e-08 | 7.279220e-07 | 7.731816e-07 | 7.473077e-07 | 6.931449e-07 | 7.097192e-07 | 5.256441e-07 | 7.511199e-07 | 7.802103e-07 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 818 | 5.753508e-12 | 1.093513e-11 | 30.640113 | 80.052673 | 1.079258e-06 | -1.289053e-06 | 5.230472e-07 | 1.035435e-06 | 3.879043e-12 | 1.241866e-11 | ... | 3.259420e-08 | 1.426668e-08 | 6.924461e-07 | 7.150275e-07 | 6.614229e-07 | 5.968342e-07 | 5.177721e-07 | 5.113420e-07 | 7.761048e-07 | 7.470384e-07 |
| 819 | 2.831619e-13 | 5.802462e-13 | 4.767611 | 26.602752 | 1.376398e-07 | -3.205892e-07 | 2.276967e-07 | 2.056310e-07 | 1.838487e-13 | 5.825916e-13 | ... | 3.056678e-08 | -7.938843e-09 | 6.200787e-07 | 8.224207e-07 | 7.316316e-07 | 6.059367e-07 | 5.789283e-07 | 4.807172e-07 | 6.892734e-07 | 6.667963e-07 |
| 820 | 6.896188e-13 | 5.856500e-13 | 32.432801 | 80.746495 | 5.683229e-09 | 3.632243e-07 | 4.151775e-07 | 1.203353e-07 | 6.570139e-13 | 2.199750e-12 | ... | 1.433694e-08 | 3.156120e-08 | 7.659014e-07 | 7.924929e-07 | 7.910302e-07 | 6.446463e-07 | 5.405371e-07 | 5.128369e-07 | 9.115879e-07 | 9.027841e-07 |
| 821 | 1.487735e-13 | 1.638185e-13 | 3.440259 | 49.056641 | -1.520831e-07 | 4.612325e-08 | 1.185922e-07 | 1.970464e-07 | 2.407891e-13 | 3.890695e-14 | ... | 2.972890e-08 | 2.393718e-08 | 6.741715e-07 | 8.067064e-07 | 7.425204e-07 | 6.147657e-07 | 5.899730e-07 | 6.015475e-07 | 8.484086e-07 | 7.875193e-07 |
| 822 | 3.172817e-12 | 2.518842e-13 | 43.820945 | 79.002367 | 7.985623e-07 | 2.192741e-07 | 3.943380e-07 | 1.220242e-07 | 3.693911e-12 | 1.030655e-13 | ... | 3.591130e-08 | 1.114778e-08 | 5.847602e-07 | 7.148684e-07 | 6.646047e-07 | 5.982145e-07 | 5.157508e-07 | 4.775573e-07 | 7.562003e-07 | 7.319647e-07 |
823 rows × 544 columns
yasa¶Detección de Espigas (Spindles). yasa proporciona herramientas para la detección de espigas, que son eventos característicos del sueño ligero. Las características extraídas de las espigas, como la densidad de espigas, la duración y la frecuencia, pueden ser útiles para identificar patrones de sueño ligero. La detección de espigas puede ayudar a identificar las etapas del sueño:
1) Distinguir entre las etapas del sueño ligero y profundo: Las espigas suelen ser más prominentes durante el sueño de ondas lentas (etapas 3 y 4 del sueño NREM), que son las etapas de sueño más profundo. La detección de espigas puede ayudar a identificar cuándo un individuo está en estas etapas de sueño más restauradoras.
2) Identificar la transición entre etapas del sueño: Las espigas pueden ocurrir durante la transición entre diferentes etapas del sueño, como el paso de la vigilia al sueño ligero o del sueño ligero al sueño profundo. La presencia de espigas en el EEG puede indicar estos cambios de etapa.
3) Monitorear la calidad del sueño: La presencia de espigas en el EEG también puede proporcionar información sobre la calidad del sueño. Una mayor cantidad de espigas puede sugerir una mayor actividad cerebral durante el sueño, lo que a menudo se asocia con una menor calidad del sueño.
4) Evaluar trastornos del sueño: La detección de espigas en el EEG también se utiliza para identificar y diagnosticar trastornos del sueño, como el síndrome de piernas inquietas o la epilepsia nocturna, que pueden afectar la calidad y la estructura del sueño.
Detección de Ondas Lentas (Slow Waves). Las ondas lentas son otro componente importante del sueño. yasa permite la detección y caracterización de estas ondas, lo que es esencial para analizar la calidad del sueño profundo. La detección de ondas lentas puede proporcionar información sobre las etapas del sueño. Permite distinguir claramente entre el sueño de ondas lentas y el sueño REM (movimiento rápido de los ojos), que son las dos principales etapas del sueño:
1) Evaluar la calidad del sueño: La cantidad y la amplitud de las ondas lentas en el EEG pueden proporcionar información sobre la calidad del sueño profundo. Un aumento en la cantidad de ondas lentas suele asociarse con un sueño de mejor calidad y un mayor grado de restauración física y mental.
2) Monitorear trastornos del sueño: La falta de ondas lentas o una alteración en su patrón puede ser indicativa de trastornos del sueño, como el insomnio o la apnea del sueño. La detección de problemas en las ondas lentas en el EEG puede ser útil en la evaluación y el diagnóstico de estos trastornos.
3) Evaluar la eficacia de tratamientos: La detección de ondas lentas puede ser útil para evaluar la eficacia de tratamientos para trastornos del sueño o para seguimiento en investigaciones relacionadas con el sueño. Los cambios en las ondas lentas pueden indicar mejoras o empeoramientos en la calidad del sueño de una persona.
Análisis de Variabilidad de la Frecuencia Cardíaca (HRV). La biblioteca también proporciona herramientas para calcular la variabilidad de la frecuencia cardíaca a partir de los datos de ECG en el PSG. Esto puede ser útil para evaluar el sistema nervioso autónomo durante el sueño. En concreto permite:
1) Distinguir entre el sueño REM y el sueño NREM: Durante el sueño, el sistema nervioso autónomo regula la frecuencia cardíaca, y esta regulación cambia de manera característica en las diferentes etapas del sueño. El análisis de la HRV puede ayudar a distinguir entre el sueño de movimiento rápido de los ojos (REM) y el sueño no REM (NREM) al evaluar patrones de variabilidad específicos en la frecuencia cardíaca.
2) Identificar la activación del sistema nervioso autónomo: Las fluctuaciones en la HRV pueden indicar la activación del sistema nervioso simpático (que tiende a acelerar la frecuencia cardíaca) o del sistema nervioso parasimpático (que tiende a desacelerarla). La HRV puede mostrar cómo estos sistemas se activan y desactivan durante el sueño, lo que es relevante para comprender las respuestas del cuerpo a los diferentes estados de sueño.
3) Evaluar la calidad del sueño: El análisis de la HRV también puede utilizarse para evaluar la calidad del sueño. Una HRV alterada o anormal puede indicar problemas en la regulación del sistema nervioso durante el sueño y puede estar asociada con trastornos del sueño, como la apnea del sueño o el insomnio.
4) Monitorizar cambios en el sueño a lo largo del tiempo: El análisis de la HRV se puede utilizar para realizar un seguimiento de las variaciones en el sueño a lo largo del tiempo. Esto es útil en investigaciones y estudios clínicos para evaluar cómo los tratamientos o intervenciones afectan la calidad del sueño y las respuestas del sistema nervioso autónomo.
(EEG-HRV coupling). El acoplamiento entre EEG y HRV puede ofrecer información sobre la interacción entre el sistema nervioso central y el sistema nervioso autónomo. Las características relacionadas con este acoplamiento podrían ser útiles para comprender la relación entre el cerebro y el corazón durante el sueño.
Identificar transiciones entre etapas del sueño: Durante la noche, las personas suelen pasar de una etapa de sueño a otra. El análisis del acoplamiento EEG-HRV puede ayudar a detectar estas transiciones al mostrar cambios en la sincronización entre la actividad cerebral y la frecuencia cardíaca. Esto puede ser útil para identificar cuándo ocurren las transiciones entre sueño ligero, sueño profundo y sueño REM.
Evaluar la calidad del sueño: La calidad del sueño se relaciona con la eficiencia de las transiciones entre las etapas del sueño y la estabilidad de estas etapas. El acoplamiento EEG-HRV puede proporcionar información sobre cómo la calidad del sueño se ve afectada por la sincronización entre la actividad cerebral y la frecuencia cardíaca.
Identificar trastornos del sueño: La presencia de patrones anómalos en el acoplamiento EEG-HRV puede ser indicativa de trastornos del sueño, como la apnea del sueño o el síndrome de piernas inquietas. Estos trastornos pueden afectar la sincronización entre la actividad cerebral y la actividad cardíaca durante el sueño.
(SO-sigma coupling). Se centra en el acoplamiento entre las oscilaciones lentas (Slow Oscillations, SO) y las oscilaciones sigma durante el sueño. El acoplamiento entre estas oscilaciones puede proporcionar información valiosa sobre la calidad del sueño y los patrones de sueño. La extracción de características relacionadas con este acoplamiento podría ser relevante para caracterizar eventos específicos durante el sueño. La detección del acoplamiento SO-sigma en el EEG puede proporcionar evidencias sobre las etapas del sueño de la siguiente manera:
1) Distinción entre el sueño REM y el sueño NREM: Las ondas sigma son prominentes durante el sueño REM, que es una etapa del sueño caracterizada por un sueño más activo y el movimiento rápido de los ojos. Las ondas delta son más prominentes en el sueño NREM, particularmente en las etapas más profundas. La presencia de acoplamiento SO-sigma puede indicar la transición entre estas dos etapas del sueño y ayudar a identificar cuándo una persona está en sueño REM.
2) Identificación de microdespertares: Durante el sueño, las personas pueden experimentar microdespertares breves que los sacan de un estado de sueño más profundo. Estos microdespertares suelen estar asociados con cambios en el acoplamiento SO-sigma, ya que se produce una transición de ondas sigma del sueño REM a ondas delta del sueño NREM.
3) Evaluación de la calidad del sueño: La presencia de acoplamiento SO-sigma anormal o irregular puede ser indicativa de problemas en la calidad del sueño. Si el patrón de acoplamiento es disruptivo o inusual, esto puede sugerir trastornos del sueño, como el insomnio o el síndrome de las piernas inquietas, que pueden afectar negativamente la calidad del sueño.
Extracción de Características de EEG. yasa permite calcular características específicas de la señal EEG, como la potencia en diferentes bandas de frecuencia (delta, theta, alpha, beta, gamma), lo que puede proporcionar información sobre la actividad cerebral durante el sueño.
1) Delta (0.5-4 Hz): Las ondas delta son de baja frecuencia y están asociadas principalmente con el sueño profundo, específicamente con las etapas 3 y 4 del sueño no REM (NREM). Un aumento en la potencia de la banda delta en el EEG indica una mayor presencia de sueño profundo, lo que sugiere una fase de sueño reparador.
2) Theta (4-8 Hz): Las ondas theta son típicas del sueño ligero y también pueden estar presentes en las etapas iniciales del sueño NREM. Un aumento en la potencia de la banda theta indica una transición desde la vigilia al sueño ligero.
3) Alpha (8-13 Hz): Las ondas alpha son comunes durante la vigilia y pueden disminuir en intensidad durante el sueño. La presencia de ondas alfa en el EEG durante el sueño generalmente se asocia con el sueño ligero.
4) Beta (13-30 Hz): Las ondas beta son típicas de la vigilia y la actividad mental consciente. Su presencia en el EEG durante el sueño sugiere actividad cerebral activa, como durante el sueño REM o sueño paradójico.
5) Gamma (30-100 Hz): Las ondas gamma son de alta frecuencia y se asocian con la activación cortical. Su presencia en el EEG puede indicar la fase de sueño REM, durante la cual se producen sueños vívidos y el cerebro muestra actividad similar a la vigilia.
(Nonlinear features). Este notebook se enfoca en la extracción de características no lineales a partir de datos de EEG. Las características no lineales pueden proporcionar información adicional sobre la complejidad de la actividad cerebral y pueden ser relevantes para la detección de patrones anómalos o para caracterizar eventos específicos en la señal EEG. Estas características no lineales pueden proporcionar información adicional sobre las etapas del sueño y otros aspectos del sueño que no se pueden capturar completamente con el análisis de características lineales. Algunas formas en que el análisis de características no lineales puede aportar a la comprensión de las etapas del sueño son las siguientes:
1) Detección de transiciones entre etapas del sueño: Las características no lineales pueden ayudar a identificar las transiciones entre las diferentes etapas del sueño, como la transición del sueño REM (movimiento rápido de los ojos) al sueño NREM (movimiento no rápido de los ojos) y viceversa. Estas transiciones pueden tener patrones no lineales específicos que se pueden detectar mediante análisis avanzados.
2) Evaluación de la variabilidad temporal: Las características no lineales pueden proporcionar información sobre la variabilidad temporal de las señales del EEG durante el sueño. Esta variabilidad temporal puede reflejar cambios en la estabilidad de las etapas del sueño y la calidad del sueño.
3) Identificación de patrones no convencionales: Las características no lineales pueden ayudar a identificar patrones no convencionales en las señales del EEG que pueden estar relacionados con condiciones específicas del sueño o trastornos del sueño. Esto es particularmente relevante en la detección de trastornos del sueño, como la apnea del sueño o el síndrome de piernas inquietas.
4) Evaluación de la coherencia y la conectividad: Las características no lineales también se utilizan para evaluar la coherencia y la conectividad entre diferentes áreas del cerebro durante el sueño. Esto puede proporcionar información sobre cómo se comunican diferentes regiones cerebrales durante las etapas del sueño.
Seguidamente, procedemos a calcular las características con una función que lo engloba todo (segmenta en ventanas de 30 segundos y calcula las características del API).
import yasa
print("Versión de YASA instalada:", yasa.__version__)
Versión de YASA instalada: 0.6.4
# Obtenemos algunas características utilizando la función SleepStaging y
for sub in pp.lSubject:
print('Subject:', sub.subject)
raw = sub.get_raw_eData()
sls = yasa.SleepStaging(raw, eeg_name = 'C4-A1',
eog_name = 'LOC-A2',
emg_name = 'X1')
# Eliminamos las 30 últimas épocas:
yasa_features = sls.get_features()[:-30]
sub.set_yasa_feature_norm_matrix(yasa_features)
pp.lSubject[0].get_yasa_feature_norm_matrix()
Subject: 6 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 4 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 8 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 3 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 7 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 5 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 2 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 9 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 1 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...). Subject: 10 NOTE: pick_channels() is a legacy function. New code should use inst.pick(...).
| eeg_abspow | eeg_abspow_c7min_norm | eeg_abspow_p2min_norm | eeg_alpha | eeg_alpha_c7min_norm | eeg_alpha_p2min_norm | eeg_at | eeg_at_c7min_norm | eeg_at_p2min_norm | eeg_beta | ... | eog_skew_c7min_norm | eog_skew_p2min_norm | eog_std | eog_std_c7min_norm | eog_std_p2min_norm | eog_theta | eog_theta_c7min_norm | eog_theta_p2min_norm | time_hour | time_norm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| epoch | |||||||||||||||||||||
| 0 | 5.399343 | 0.148168 | 0.907698 | 0.104052 | 1.273231 | -0.145396 | 1.747386 | 1.809266 | 1.496076 | 0.137878 | ... | -0.231166 | -0.057000 | 4.256332 | 0.428076 | 2.204360 | 0.014738 | 0.108276 | -0.771056 | 0.000000 | 0.000000 |
| 1 | 1.546661 | 0.079136 | 0.459812 | 0.488675 | 1.312799 | 0.871101 | 1.988468 | 1.738572 | 1.665018 | 0.082483 | ... | -0.209726 | -0.121376 | 0.972484 | 0.257628 | 1.052684 | 0.270484 | 0.150062 | -0.024735 | 0.008333 | 0.001174 |
| 2 | 1.463680 | 0.027225 | 0.304086 | 0.415958 | 1.290048 | 1.081813 | 1.551462 | 1.656668 | 1.517173 | 0.098716 | ... | -0.162498 | -0.219808 | 0.903806 | 0.158585 | 0.652734 | 0.229639 | 0.138346 | 0.144575 | 0.016667 | 0.002347 |
| 3 | 1.002520 | -0.015675 | 0.199417 | 0.392485 | 1.266602 | 1.156152 | 1.756101 | 1.605253 | 1.514953 | 0.090743 | ... | -0.101075 | -0.275230 | 1.090249 | 0.084518 | 0.485453 | 0.167393 | 0.115337 | 0.138408 | 0.025000 | 0.003521 |
| 4 | 1.002217 | -0.051649 | -0.056173 | 0.456902 | 1.218936 | 1.622415 | 1.982643 | 1.528801 | 1.597383 | 0.090143 | ... | -0.018408 | -0.174047 | 0.836021 | 0.020465 | -0.114314 | 0.181707 | 0.092402 | 0.382034 | 0.033333 | 0.004695 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 818 | 1.028830 | -0.108890 | -0.080878 | 0.209690 | 0.597308 | 0.519031 | 0.967755 | 0.665006 | 0.537351 | 0.026888 | ... | -0.373300 | -0.244966 | 0.822869 | -0.285017 | -0.181453 | 0.146267 | 0.099286 | 0.130593 | 6.816667 | 0.960094 |
| 819 | 1.052436 | -0.103241 | -0.104524 | 0.212754 | 0.599079 | 0.568998 | 1.109394 | 0.684633 | 0.578471 | 0.034760 | ... | -0.471817 | -0.187725 | 0.798703 | -0.281959 | -0.196245 | 0.171105 | 0.093444 | 0.175433 | 6.825000 | 0.961268 |
| 820 | 1.359540 | -0.097897 | -0.085965 | 0.313967 | 0.606365 | 0.577785 | 1.807835 | 0.699035 | 0.813016 | 0.026110 | ... | -0.551604 | -0.229835 | 0.801294 | -0.280122 | -0.201357 | 0.124474 | 0.085002 | -0.008229 | 6.833333 | 0.962441 |
| 821 | 1.063286 | -0.094848 | -0.085875 | 0.211460 | 0.591722 | 0.557154 | 0.736069 | 0.664120 | 0.666199 | 0.036550 | ... | -0.583051 | -0.441468 | 0.959299 | -0.276013 | -0.188032 | 0.165690 | 0.092969 | 0.029388 | 6.841667 | 0.963615 |
| 822 | 1.012121 | -0.087669 | -0.086846 | 0.219702 | 0.579714 | 0.570384 | 0.904733 | 0.647850 | 0.644117 | 0.038175 | ... | -0.546312 | -0.403118 | 0.711851 | -0.259666 | -0.207500 | 0.195659 | 0.102436 | 0.101456 | 6.850000 | 0.964789 |
823 rows × 149 columns
| abspow | abspow_c7min_norm | abspow_p2min_norm | alpha | alpha_c7min_norm | alpha_p2min_norm | at | at_c7min_norm | at_p2min_norm | beta | ... | skew_c7min_norm | skew_p2min_norm | std | std_c7min_norm | std_p2min_norm | theta | theta_c7min_norm | theta_p2min_norm | time_hour | time_norm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| eeg | 2.217383 | 0.207680 | 0.001994 | 0.306759 | 0.378277 | 0.850549 | 4.186395 | 0.938804 | 2.029416 | 0.126715 | ... | -0.450356 | -0.518535 | 3.657893 | 1.719164 | 2.064860 | 0.026997 | -0.751876 | -0.772705 | 0.000000 | 0.000000 |
| eog | 2.517489 | 0.033908 | 0.055935 | 0.134619 | 0.103857 | 0.145556 | 1.014488 | 0.047966 | 0.122808 | 0.021424 | ... | -0.127470 | -0.462168 | 0.920967 | -0.285169 | -0.200311 | 0.166790 | 0.272841 | 0.290148 | 7.658333 | 0.964323 |
| emg | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
Normalizamos las características definidas en el artículo de referencia [1].
# Comprobamos si la matriz tiene valores Nan, inf o nulos. En ese caso se le asignaria el valor 0:
# valores_nan = paper_features_0.isna().sum().sum()
# Mostrar la cantidad de valores NaN
# print("Cantidad de valores NaN: ", valores_nan)
for sub in pp.lSubject:
print ('Subject:', sub.subject)
paper_features_norm = extract.normalice_matrix_features(sub.get_paper_feature_norm_matrix())
sub.set_paper_feature_norm_matrix(paper_features_norm)
Subject: 6 Subject: 4 Subject: 8 Subject: 3 Subject: 7 Subject: 5 Subject: 2 Subject: 9 Subject: 1 Subject: 10
yasa¶Normalizamos las características creadas por la función yasa (yasa.SleepStaging(raw, eeg_name = 'C4-A1', eog_name = 'LOC-A2', emg_name = 'X1')).
# Comprobamos si la matriz tiene valores Nan, inf o nulos. En ese caso se le asignaria el valor 0.
# valores_nan = yasa_features.isna().sum().sum()
# Mostrar la cantidad de valores NaN
# print("Cantidad de valores NaN: ", valores_nan)
for sub in pp.lSubject:
print('Subject:', sub.subject)
yasa_features_norm = extract.normalice_matrix_features(sub.get_yasa_feature_norm_matrix())
sub.set_yasa_feature_norm_matrix(yasa_features_norm)
Subject: 6 Subject: 4 Subject: 8 Subject: 3 Subject: 7 Subject: 5 Subject: 2 Subject: 9 Subject: 1 Subject: 10
Dividimos los sujetos en entrenamiento, validación y prueba (6, 2, 2). Del conjunto de 10 sujetos, tendremos 6 sujetos para entrenamiento, 2 para validación y otros 2 para prueba.
X_train_o, X_val_o, X_test_o = pp.split_data(0.6, 0.2, 0.2)
print(' Sujetos para Entrenamiento:')
for i in range(len(X_train_o)):
print('\t\t Subject:', X_train_o[i].subject)
print(' Sujetos para validación:')
for i in range(len(X_val_o)):
print ('\t\t Subject:', X_val_o[i].subject)
print(' Sujetos para test:')
for i in range(len(X_test_o)):
print('\t\t Subject:', X_test_o[i].subject)
Sujetos para Entrenamiento: Subject: 5 Subject: 4 Subject: 6 Subject: 3 Subject: 7 Subject: 1 Sujetos para validación: Subject: 8 Subject: 10 Sujetos para test: Subject: 9 Subject: 2
Para realizar una selección de características, nos basamos en lo predicho por los expertos. De esta forma obtenemos una primera selección óptima con la prueba f-anova y cross-validation que se va optimizando iterativamente eliminando aquella característica que aporta menos al conjunto obtenido en primera insancia. Así hasta encontrar un peor valor, en este momento ya tenemos el listado de características óptimas para alcanzar lo predicho por los expertos.
for sub in pp.lSubject:
print ('Subject:',sub.subject)
# Cargamos las etiquetas del primer médico y eliminamos las 30 por que lo indica el artículo
path_lab = f'data/{sub.subject}/{sub.subject}_1.txt' # Camino a las etiquetas de marcado del clínico 1
labels1 = pd.read_csv(path_lab, header = None).squeeze("columns")
# Atención con la codificación si se quiere utilizar las funciones de YASA...
labels1[labels1 == 5] = 4
labels1 = labels1[:-30]
print(f'labels1:{len (labels1)}')
# Cargamos las etiquetas del primer médico y eliminamos las 30 por que lo indica el artículo
path_lab = f'data/{sub.subject}/{sub.subject}_2.txt' # Camino a las etiquetas de marcado del clínico 1
labels2 = pd.read_csv(path_lab, header = None).squeeze("columns")
# Atención con la codificación si se quiere utilizar las funciones de YASA...
labels2 [labels2 == 5] = 4
labels2 = labels2[:-30]
print (f'labels2:{len(labels2)}')
differences = np.sum(np.array(labels2) == np.array(labels1))
if differences != len(labels2) :
print("Las listas tienen elementos diferentes.")
else:
print("Las listas son iguales.")
sub.set_labels ([labels1, labels2])
Subject: 5 labels1:914 labels2:914 Las listas tienen elementos diferentes. Subject: 4 labels1:764 labels2:764 Las listas tienen elementos diferentes. Subject: 6 labels1:823 labels2:823 Las listas tienen elementos diferentes. Subject: 3 labels1:794 labels2:794 Las listas tienen elementos diferentes. Subject: 7 labels1:784 labels2:784 Las listas tienen elementos diferentes. Subject: 1 labels1:924 labels2:924 Las listas tienen elementos diferentes. Subject: 8 labels1:970 labels2:970 Las listas son iguales. Subject: 10 labels1:766 labels2:766 Las listas son iguales. Subject: 9 labels1:939 labels2:939 Las listas son iguales. Subject: 2 labels1:911 labels2:911 Las listas tienen elementos diferentes.
l_9_0 = pp.lSubject[9].get_labels() [0].to_list()
l_9_1 = pp.lSubject[9].get_labels() [1].to_list()
l_4_0 = pp.lSubject[4].get_labels() [0].to_list()
l_4_1 = pp.lSubject[4].get_labels() [1].to_list()
# Verificamos si hay elementos diferentes entre las dos listas
print(len(l_9_0))
print(len(l_4_0))
print(np.sum(np.array(l_4_0) == np.array(l_4_1)))
differences = np.sum(np.array(l_4_1) == np.array(l_4_0))
if differences != len(l_4_0) :
print("Las listas tienen elementos diferentes.")
else:
print("Las listas son iguales.")
911 784 642 Las listas tienen elementos diferentes.
Para la selección de características de forma automática utilizamos los datos X_train_o e y_train.
tupla_paper_features_labels = []
tupla_yasa_features_labels = []
for sub in X_train_o:
tupla_paper_features_labels.extend(sub.get_tupla_paper_feeatures_labels())
tupla_yasa_features_labels.extend(sub.get_tupla_yasa_feeatures_labels())
Solo los sujetos dentro del conjunto de entrenamiento:
selected_paper_features_for_all_subjects = extract.select_features_by_consensus(tupla_paper_features_labels, threshold = 8)
_ = [sub.set_paper_optim_features(selected_paper_features_for_all_subjects) for sub in pp.lSubject]
print(len(selected_paper_features_for_all_subjects))
print(len(pp.lSubject[4].get_paper_feature_norm_matrix().columns))
print(set(range(0, len(pp.lSubject[4].get_paper_feature_norm_matrix().columns))) - set(selected_paper_features_for_all_subjects))
410
544
{512, 513, 2, 514, 515, 516, 517, 518, 519, 9, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 409, 411, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511}
pp.lSubject[4].get_paper_optim_features()
| Energy_s_nivel-1_Canal-O1-A2 | Energy_d_nivel-1_Canal-O1-A2 | Percentage_energy_d_nivel-1_Canal-O1-A2 | Mean_s_nivel-1_Canal-O1-A2 | Mean_d_nivel-1_Canal-O1-A2 | Std_deviation_s_nivel-1_Canal-O1-A2 | Std_deviation_d_nivel-1_Canal-O1-A2 | Energy_s_nivel-1_Canal-F4-A1 | Percentage_energy_s_nivel-1_Canal-F4-A1 | Percentage_energy_d_nivel-1_Canal-F4-A1 | ... | fc_Beta 1_Canal-F3-A2 | fr_Beta 1_Canal-F3-A2 | sfc_Beta 1_Canal-F3-A2 | fc_Beta 1_Canal-C3-A2 | fr_Beta 1_Canal-C3-A2 | sfc_Beta 1_Canal-C3-A2 | fc_Beta 2_Canal-O1-A2 | sfc_Beta 2_Canal-O1-A2 | fr_Beta 2_Canal-F4-A1 | sfc_Beta 2_Canal-F4-A1 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.121817 | 0.164942 | 1.0 | 0.000000 | 0.000000 | 0.328222 | 0.394426 | 0.118549 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 1 | 0.089184 | 0.000000 | 1.0 | 0.000000 | 0.000000 | 0.307531 | 0.336302 | 0.106424 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 2 | 0.000000 | 0.000000 | 1.0 | 0.179112 | 0.333611 | 0.182983 | 0.340291 | 0.000000 | 0.387464 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 3 | 0.053080 | 0.000000 | 1.0 | 0.000000 | 0.299163 | 0.252827 | 0.339209 | 0.050966 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 4 | 0.000000 | 0.105655 | 1.0 | 0.000000 | 0.418551 | 0.168176 | 0.388921 | 0.000000 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 779 | 0.146112 | 0.000000 | 1.0 | 0.359264 | 0.000000 | 0.326687 | 0.359162 | 0.144081 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.430525 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 780 | 0.143906 | 0.000000 | 1.0 | 0.198947 | 0.000000 | 0.385491 | 0.390016 | 0.149802 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 781 | 0.102908 | 0.000000 | 1.0 | 0.000000 | 0.000000 | 0.271672 | 0.329568 | 0.061243 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 782 | 0.101032 | 0.000000 | 1.0 | 0.253172 | 0.000000 | 0.334455 | 0.257739 | 0.096737 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 783 | 0.035676 | 0.000000 | 1.0 | 0.000000 | 0.000000 | 0.204819 | 0.334520 | 0.035162 | 1.000000 | 0.0 | ... | 0.0 | 0.0 | 0.000000 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
784 rows × 410 columns
print (pp.lSubject [4].get_paper_optim_features ().columns)
features_rejected = set (range(0,len(pp.lSubject [4].get_paper_feature_norm_matrix ().columns)))-set(selected_paper_features_for_all_subjects)
print (features_rejected)
features_rejected = list (features_rejected)
for i in range (len(features_rejected)):
print (pp.lSubject [4].get_paper_feature_norm_matrix ().columns [features_rejected [i]])
Index(['Energy_s_nivel-1_Canal-O1-A2', 'Energy_d_nivel-1_Canal-O1-A2',
'Percentage_energy_d_nivel-1_Canal-O1-A2', 'Mean_s_nivel-1_Canal-O1-A2',
'Mean_d_nivel-1_Canal-O1-A2', 'Std_deviation_s_nivel-1_Canal-O1-A2',
'Std_deviation_d_nivel-1_Canal-O1-A2', 'Energy_s_nivel-1_Canal-F4-A1',
'Percentage_energy_s_nivel-1_Canal-F4-A1',
'Percentage_energy_d_nivel-1_Canal-F4-A1',
...
'fc_Beta 1_Canal-F3-A2', 'fr_Beta 1_Canal-F3-A2',
'sfc_Beta 1_Canal-F3-A2', 'fc_Beta 1_Canal-C3-A2',
'fr_Beta 1_Canal-C3-A2', 'sfc_Beta 1_Canal-C3-A2',
'fc_Beta 2_Canal-O1-A2', 'sfc_Beta 2_Canal-O1-A2',
'fr_Beta 2_Canal-F4-A1', 'sfc_Beta 2_Canal-F4-A1'],
dtype='object', length=410)
{512, 513, 2, 514, 515, 516, 517, 518, 519, 9, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 409, 411, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511}
Std_AR_Canal_O1-A2
Std_AR_Canal_F4-A1
Percentage_energy_s_nivel-1_Canal-O1-A2
Std_AR_Canal_C4-A1
Std_AR_Canal_O2-A1
Std_AR_Canal_ROC-A1
Std_AR_Canal_LOC-A2
Std_AR_Canal_F3-A2
Std_AR_Canal_C3-A2
Energy_d_nivel-1_Canal-F4-A1
Percentil25_O1-A2
Percentil25_F4-A1
Percentil25_C4-A1
Percentil25_O2-A1
Percentil25_ROC-A1
Percentil25_LOC-A2
Percentil25_F3-A2
Percentil25_C3-A2
Percentil50_O1-A2
Percentil50_F4-A1
Percentil50_C4-A1
Percentil50_O2-A1
Percentil50_ROC-A1
Percentil50_LOC-A2
Percentil50_F3-A2
Percentil50_C3-A2
Percentil75_O1-A2
Percentil75_F4-A1
Percentil75_C4-A1
Percentil75_O2-A1
Percentil75_ROC-A1
Percentil75_LOC-A2
Percentil75_F3-A2
Percentil75_C3-A2
fr_Beta 2_Canal-O1-A2
fc_Beta 2_Canal-F4-A1
fc_Beta 2_Canal-C4-A1
fr_Beta 2_Canal-C4-A1
sfc_Beta 2_Canal-C4-A1
fc_Beta 2_Canal-O2-A1
fr_Beta 2_Canal-O2-A1
sfc_Beta 2_Canal-O2-A1
fc_Beta 2_Canal-ROC-A1
fr_Beta 2_Canal-ROC-A1
sfc_Beta 2_Canal-ROC-A1
fc_Beta 2_Canal-LOC-A2
fr_Beta 2_Canal-LOC-A2
sfc_Beta 2_Canal-LOC-A2
fc_Beta 2_Canal-F3-A2
fr_Beta 2_Canal-F3-A2
sfc_Beta 2_Canal-F3-A2
fc_Beta 2_Canal-C3-A2
fr_Beta 2_Canal-C3-A2
sfc_Beta 2_Canal-C3-A2
Activity_Canal_O1-A2
Activity_Canal_F4-A1
Activity_Canal_C4-A1
Activity_Canal_O2-A1
Activity_Canal_ROC-A1
Activity_Canal_LOC-A2
Activity_Canal_F3-A2
Activity_Canal_C3-A2
Mobility_Canal_O1-A2
Mobility_Canal_F4-A1
Mobility_Canal_C4-A1
Mobility_Canal_O2-A1
Mobility_Canal_ROC-A1
Mobility_Canal_LOC-A2
Mobility_Canal_F3-A2
Mobility_Canal_C3-A2
Complexity_Canal_O1-A2
Complexity_Canal_F4-A1
Complexity_Canal_C4-A1
Complexity_Canal_O2-A1
Complexity_Canal_ROC-A1
Complexity_Canal_LOC-A2
Complexity_Canal_F3-A2
Complexity_Canal_C3-A2
Skewness_Canal_O1-A2
Skewness_Canal_F4-A1
Skewness_Canal_C4-A1
Skewness_Canal_O2-A1
Skewness_Canal_ROC-A1
Skewness_Canal_LOC-A2
Skewness_Canal_F3-A2
Skewness_Canal_C3-A2
Kurtosis_Canal_O1-A2
Kurtosis_Canal_F4-A1
Kurtosis_Canal_C4-A1
Kurtosis_Canal_O2-A1
Kurtosis_Canal_ROC-A1
Kurtosis_Canal_LOC-A2
Kurtosis_Canal_F3-A2
Kurtosis_Canal_C3-A2
AR_Coef_0_Canal_O1-A2
AR_Coef_0_Canal_F4-A1
AR_Coef_0_Canal_C4-A1
AR_Coef_0_Canal_O2-A1
AR_Coef_0_Canal_ROC-A1
AR_Coef_0_Canal_LOC-A2
AR_Coef_0_Canal_F3-A2
AR_Coef_0_Canal_C3-A2
AR_Coef_1_Canal_O1-A2
AR_Coef_1_Canal_F4-A1
AR_Coef_1_Canal_C4-A1
AR_Coef_1_Canal_O2-A1
AR_Coef_1_Canal_ROC-A1
AR_Coef_1_Canal_LOC-A2
AR_Coef_1_Canal_F3-A2
AR_Coef_1_Canal_C3-A2
AR_Coef_2_Canal_O1-A2
AR_Coef_2_Canal_F4-A1
AR_Coef_2_Canal_C4-A1
AR_Coef_2_Canal_O2-A1
AR_Coef_2_Canal_ROC-A1
AR_Coef_2_Canal_LOC-A2
AR_Coef_2_Canal_F3-A2
AR_Coef_2_Canal_C3-A2
AR_Coef_3_Canal_O1-A2
AR_Coef_3_Canal_F4-A1
AR_Coef_3_Canal_C4-A1
AR_Coef_3_Canal_O2-A1
AR_Coef_3_Canal_ROC-A1
AR_Coef_3_Canal_LOC-A2
AR_Coef_3_Canal_F3-A2
AR_Coef_3_Canal_C3-A2
Media_AR_Canal_O1-A2
Media_AR_Canal_F4-A1
Media_AR_Canal_C4-A1
Media_AR_Canal_O2-A1
Media_AR_Canal_ROC-A1
Media_AR_Canal_LOC-A2
Media_AR_Canal_F3-A2
Media_AR_Canal_C3-A2
yasaSólo los sujetos dentro del grupo de entrenamiento:
selected_yasa_features_for_all_subjects = extract.select_features_by_consensus(tupla_yasa_features_labels, threshold = 8)
_ = [sub.set_yasa_optim_features(selected_yasa_features_for_all_subjects) for sub in pp.lSubject]
print(len(selected_yasa_features_for_all_subjects))
print(len(pp.lSubject[4].get_yasa_feature_norm_matrix().columns))
print(set(range(0, len(pp.lSubject[4].get_yasa_feature_norm_matrix().columns))) - set(selected_yasa_features_for_all_subjects))
133
149
{135, 136, 137, 138, 139, 12, 140, 141, 142, 143, 144, 145, 146, 147, 148, 29}
print (pp.lSubject [4].get_paper_optim_features ().columns)
features_rejected = set (range(0,len(pp.lSubject [4].get_yasa_feature_norm_matrix ().columns)))-set(selected_yasa_features_for_all_subjects)
print (features_rejected)
features_rejected = list (features_rejected)
for i in range (len(features_rejected)):
print (pp.lSubject [4].get_yasa_feature_norm_matrix ().columns [features_rejected [i]])
Index(['Energy_s_nivel-1_Canal-O1-A2', 'Energy_d_nivel-1_Canal-O1-A2',
'Percentage_energy_d_nivel-1_Canal-O1-A2', 'Mean_s_nivel-1_Canal-O1-A2',
'Mean_d_nivel-1_Canal-O1-A2', 'Std_deviation_s_nivel-1_Canal-O1-A2',
'Std_deviation_d_nivel-1_Canal-O1-A2', 'Energy_s_nivel-1_Canal-F4-A1',
'Percentage_energy_s_nivel-1_Canal-F4-A1',
'Percentage_energy_d_nivel-1_Canal-F4-A1',
...
'fc_Beta 1_Canal-F3-A2', 'fr_Beta 1_Canal-F3-A2',
'sfc_Beta 1_Canal-F3-A2', 'fc_Beta 1_Canal-C3-A2',
'fr_Beta 1_Canal-C3-A2', 'sfc_Beta 1_Canal-C3-A2',
'fc_Beta 2_Canal-O1-A2', 'sfc_Beta 2_Canal-O1-A2',
'fr_Beta 2_Canal-F4-A1', 'sfc_Beta 2_Canal-F4-A1'],
dtype='object', length=410)
{135, 136, 137, 138, 139, 12, 140, 141, 142, 143, 144, 145, 146, 147, 148, 29}
eog_sigma
eog_sigma_c7min_norm
eog_sigma_p2min_norm
eog_skew
eog_skew_c7min_norm
eeg_db
eog_skew_p2min_norm
eog_std
eog_std_c7min_norm
eog_std_p2min_norm
eog_theta
eog_theta_c7min_norm
eog_theta_p2min_norm
time_hour
time_norm
eeg_higuchi_p2min_norm
Creamos el modelo RandomForestClassifier y lo entrenamos con los 6 sujetos escogidos para entrenamiento. Obtenemos la curva ROC para entrenamiento.
from sklearn.preprocessing import label_binarize
def binarize_y (vec):
return label_binarize(vec, classes = [0, 1, 2, 3, 4])
X_train = [sub.get_paper_optim_features().values for sub in X_train_o]
X_val = [sub.get_paper_optim_features().values for sub in X_val_o]
X_test = [sub.get_paper_optim_features().values for sub in X_test_o]
y_train_0 = [binarize_y(sub.get_labels()[0]) for sub in X_train_o]
y_train_1 = [binarize_y(sub.get_labels()[1]) for sub in X_train_o]
y_val_0 = [binarize_y(sub.get_labels()[0]) for sub in X_val_o]
y_val_1 = [binarize_y(sub.get_labels()[1]) for sub in X_val_o]
y_test_0 = [binarize_y(sub.get_labels()[0]) for sub in X_test_o]
y_test_1 = [binarize_y(sub.get_labels()[1]) for sub in X_test_o]
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import label_binarize
from sklearn.multiclass import OneVsRestClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_curve, auc
# Combinamos las características y etiquetas de los seis sujetos en una sola estructura
combined_X = []
combined_y_0 = []
combined_y_1 = []
for i in range(6):
combined_X.extend(X_train[i])
combined_y_0.extend(y_train_0[i])
combined_y_1.extend(y_train_1[i])
combined_X = np.array(combined_X)
combined_y_0 = np.array(combined_y_0)
combined_y_1 = np.array(combined_y_1)
combined_x_all = np.concatenate([combined_X, combined_X])
combined_y_all = np.concatenate([combined_y_0, combined_y_1])
# Definimos el modelo de Random Forest
classifier = OneVsRestClassifier(
RandomForestClassifier(n_estimators = 1000,
criterion = "gini",
random_state = 0)
)
classifier.fit(combined_x_all, combined_y_all)
OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))RandomForestClassifier(n_estimators=1000, random_state=0)
RandomForestClassifier(n_estimators=1000, random_state=0)
from sklearn.inspection import permutation_importance
for i in range (len(X_test)):
for j in range (len(y_test_0)):
# Calcular la importancia de la permutación
result = permutation_importance(classifier, X_test [i], y_test_0 [j], n_repeats = 30, random_state = 42)
# Obtener las importancias
importances = result.importances_mean
# Obtener índices ordenados por importancia
indices = importances.argsort()[::-1]
print (f'Las 30 features más importantes para el subject {i} del Experto 1:')
# Imprimir las características ordenadas por importancia
for f in range(X_test[i].shape [1]):
print(f"\t\t{X_test[i].columns[indices[f]]}: {importances[indices[f]]}")
# Como tarda bastante, la traza que saccamos es bajo el subject 0 y experto 1 . no se hace el for...
# si se desea para los dos subjects hay que esperar bastante y consumir muchisimos recursos ciputacionales...
from sklearn.inspection import permutation_importance
for i in range (len(X_test)):
for j in range (len(y_test_0)):
# Calcular la importancia de la permutación
result = permutation_importance(classifier, X_test [i], y_test_1 [j], n_repeats = 30, random_state = 42)
# Obtener las importancias
importances = result.importances_mean
# Obtener índices ordenados por importancia
indices = importances.argsort()[::-1]
print (f' Las 30 features más importantes para el subject {i} del Experto 2:)
# Imprimir las características ordenadas por importancia
for f in range(X_test[i] .shape[1]):
print(f"{X_test[i].columns[indices[f]]}: {importances[indices[f]]}")
import numpy as np
def ajustar_umbral(probabilidades, umbral):
return [1 if prob >= umbral else 0 for prob in probabilidades]
# Ejemplo de uso
# Obtén las probabilidades para las clases en el conjunto de prueba
for i in range (len(X_test)):
probabilidades_clases = classifier.predict_proba(X_test [i])
print (probabilidades_clases)
umbral_ajustado = 0.6
predicciones_ajustadas = ajustar_umbral(probabilidades_clases, umbral_ajustado)
print(f' Las predicciones ajustadas para el subject {i}')
print(predicciones_ajustadas)
from sklearn.metrics import precision_recall_curve, auc
import matplotlib.pyplot as plt
def ajustar_umbral(probabilidades, y_true, umbral_inicial=0.5):
precision, recall, umbrales = precision_recall_curve(y_true, probabilidades)
area_bajo_curva = auc(recall, precision)
mejor_umbral = umbral_inicial
mejor_f1 = 0
for i in range(len(umbrales)):
umbral_actual = umbrales[i]
y_pred_ajustado = (probabilidades >= umbral_actual).astype(int)
f1_actual = 2 * (precision[i] * recall[i]) / (precision[i] + recall[i])
if f1_actual > mejor_f1:
mejor_f1 = f1_actual
mejor_umbral = umbral_actual
return mejor_umbral
# Ejemplo de uso
probabilidades_clases = classifier.predict_proba(X_test [0])
print (probabilidades_clases)
umbral_optimo = ajustar_umbral(probabilidades_clases, y_test_0 [0])
print("Umbral óptimo para Experto 1 del subject 1:", umbral_optimo)
umbral_optimo = ajustar_umbral(probabilidades_clases, y_test_1 [0])
print("Umbral óptimo para Experto 2 del subject 1:", umbral_optimo)
Notar que usamos directamente las probabiliades que nos proporciona el modelo. Esto es así porque las predicciones y los Auc de las curvas Roc tanto para validación como para Test son bastantes aceptables. Además, las matrices de confusión son igualmente representativas y no es necesario ajustar mucho los umbrales de probabilidad para poder predecir mejor. Si acaso, para la clase N2 podría ser que bajado algo el umbral se predijera y seleciconará mejor dicha clase.
Entrenamiento del modelo.
pp.calc_roc_acu(classifier, combined_x_all, combined_y_all)
{0: 0.9992726157134799,
1: 0.9881299416149647,
2: 0.9924173893957493,
3: 0.9981082199322366,
4: 0.9985678162262595}
Matriz de confusión para el conjunto de entrenamiento (combinado).
cm_entrenamiento = pp.calculate_multilabel_confusion_matrix(classifier, combined_x_all, combined_y_all)
print("Matriz de Confusión para el Conjunto de Entrenamiento:")
print(cm_entrenamiento)
Matriz de Confusión para el Conjunto de Entrenamiento: [[[8118 75] [ 72 1741]] [[8882 46] [ 432 646]] [[6426 287] [ 292 3001]] [[7338 95] [ 174 2399]] [[8716 41] [ 136 1113]]]
Validación.
# Realizamos la predicción para el primer sujeto de validación en X_val[0]
# Evaluamos el rendimiento en el primer sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu(classifier, X_val[0], y_val_0[0])
roc_auc_val_1 = pp.calc_roc_acu(classifier, X_val[0], y_val_1[0])
print("Exactitud para el primer sujeto y primer Experto de validación:", roc_auc_val_0)
print("Exactitud para el primer sujeto y segundo Experto de validación:", roc_auc_val_1)
print ("")
# Realizamos la predicción para el segundo sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu(classifier, X_val[0], y_val_0[0])
roc_auc_val_1 = pp.calc_roc_acu(classifier, X_val[0], y_val_1[0])
print("Exactitud para el segundo sujeto y primer Exèrto de validación:", roc_auc_val_0)
print("Exactitud para el segundo sujeto y segundo Exèrto de de validación:", roc_auc_val_1)
Exactitud para el primer sujeto y primer Experto de validación: {0: 0.994016806722689, 1: 0.8629500580720093, 2: 0.9623299500478266, 3: 0.9970827238058193, 4: 0.9867651987672588}
Exactitud para el primer sujeto y segundo Experto de validación: {0: 0.994016806722689, 1: 0.8629500580720093, 2: 0.9623299500478266, 3: 0.9970827238058193, 4: 0.9867651987672588}
Exactitud para el segundo sujeto y primer Exèrto de validación: {0: 0.994016806722689, 1: 0.8629500580720093, 2: 0.9623299500478266, 3: 0.9970827238058193, 4: 0.9867651987672588}
Exactitud para el segundo sujeto y segundo Exèrto de de validación: {0: 0.994016806722689, 1: 0.8629500580720093, 2: 0.9623299500478266, 3: 0.9970827238058193, 4: 0.9867651987672588}
Matriz de confusión para los datos de validación.
for i in range(len(X_val)):
# roc_auc_val_0 = pp.calc_roc_acu(clasificador, X_val[i], y_val_0[i])
# roc_auc_val_1 = pp.calc_roc_acu(clasificador, X_val[i], y_val_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Validación - Sujeto {i + 1}:")
cm_val_0 = pp.calculate_multilabel_confusion_matrix(classifier, X_val[i], y_val_0[i])
cm_val_1 = pp.calculate_multilabel_confusion_matrix(classifier, X_val[i], y_val_1[i])
print(f"Experto 1:\n{cm_val_0}")
print(f"Experto 2:\n{cm_val_1}")
Matriz de Confusión para el Conjunto de Validación - Sujeto 1: Experto 1: [[[581 14] [ 12 363]] [[861 0] [109 0]] [[747 29] [ 45 149]] [[816 11] [ 9 134]] [[810 11] [ 29 120]]] Experto 2: [[[581 14] [ 12 363]] [[861 0] [109 0]] [[747 29] [ 45 149]] [[816 11] [ 9 134]] [[810 11] [ 29 120]]] Matriz de Confusión para el Conjunto de Validación - Sujeto 2: Experto 1: [[[602 37] [ 5 122]] [[547 0] [217 2]] [[570 23] [100 73]] [[645 9] [ 17 95]] [[589 42] [ 37 98]]] Experto 2: [[[602 37] [ 5 122]] [[547 0] [217 2]] [[570 23] [100 73]] [[645 9] [ 17 95]] [[589 42] [ 37 98]]]
Conjunto de prueba.
# Evaluamos el rendimiento en el primer sujeto de test
roc_auc_test_0 = pp.calc_roc_acu(classifier, X_test[0], y_test_0[0])
roc_auc_test_1 = pp.calc_roc_acu(classifier, X_test[0], y_test_1[0])
print("Exactitud para el primer sujeto y primer Experto de Test:", roc_auc_test_0)
print("Exactitud para el primer sujeto y segundo Experto de Test:", roc_auc_test_1)
print ("")
# Realizamos la predicción para el segundo sujeto de test
roc_auc_test_0 = pp.calc_roc_acu(classifier, X_test[1], y_test_0[1])
roc_auc_test_1 = pp.calc_roc_acu(classifier, X_test[1], y_test_1[1])
print("Exactitud para el segundo sujeto y primer Experto de test:", roc_auc_test_0)
print("Exactitud para el segundo sujeto y segundo Experto de test:", roc_auc_test_1)
Exactitud para el primer sujeto y primer Experto de Test: {0: 0.9802054698316511, 1: 0.8291241351493904, 2: 0.8936152896004289, 3: 0.9837659508247745, 4: 0.9483764586504313}
Exactitud para el primer sujeto y segundo Experto de Test: {0: 0.9802054698316511, 1: 0.8291241351493904, 2: 0.8936152896004289, 3: 0.9837659508247745, 4: 0.9483764586504313}
Exactitud para el segundo sujeto y primer Experto de test: {0: 0.9895067817509249, 1: 0.8630192502532927, 2: 0.9134601208904614, 3: 0.9695929133074549, 4: 0.9792904073587385}
Exactitud para el segundo sujeto y segundo Experto de test: {0: 0.9951940366822506, 1: 0.8898277276456112, 2: 0.9102647924158223, 3: 0.9365774623592322, 4: 0.9742450582705556}
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [0])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_0 [0].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 9 - Expert 1')
plt.legend(loc='lower right')
plt.show()
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [1])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_1 [1].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 2 - Expert 2')
plt.legend(loc='lower right')
plt.show()
Matriz de confusión para los datos de prueba.
for i in range(len(X_test)):
# roc_auc_test_0 = pp.calc_roc_acu(classifier_0, X_test[i], y_test_0[i])
# roc_auc_test_1 = pp.calc_roc_acu(classifier_1, X_test[i], y_test_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Prueba - Sujeto {i + 1}:")
cm_test_0 = pp.calculate_multilabel_confusion_matrix(classifier, X_test[i], y_test_0[i])
cm_test_1 = pp.calculate_multilabel_confusion_matrix(classifier, X_test[i], y_test_1[i])
print(f"Experto 1:\n{cm_test_0}")
print(f"Experto 2:\n{cm_test_1}")
Matriz de Confusión para el Conjunto de Prueba - Sujeto 1: Experto 1: [[[799 18] [ 31 91]] [[771 1] [164 3]] [[495 82] [107 255]] [[654 60] [ 15 210]] [[857 19] [ 27 36]]] Experto 2: [[[799 18] [ 31 91]] [[771 1] [164 3]] [[495 82] [107 255]] [[654 60] [ 15 210]] [[857 19] [ 27 36]]] Matriz de Confusión para el Conjunto de Prueba - Sujeto 2: Experto 1: [[[803 8] [ 15 85]] [[769 1] [135 6]] [[539 49] [108 215]] [[678 36] [ 37 160]] [[755 6] [ 48 102]]] Experto 2: [[[816 22] [ 2 71]] [[803 2] [101 5]] [[535 54] [112 210]] [[645 23] [ 70 173]] [[740 4] [ 63 104]]]
import matplotlib.pyplot as plt
import numpy as np
# Obtener señales EEG para un sujeto específico (cambia el índice según se desee)
# Notar que en la ejecución teniamos el subject 9 y 2 como subjects de test
eeg_signals = pp.subject[9].getEEG()
# Obtener probabilidades de clase del clasificador
proba_predictions = classifier.predict_proba(X_test) # Ajusta X_test según tus datos de prueba
# Definir umbrales para la clasificación de las cinco clases
thresholds = [0.5, 0.2, 0.4, 0.5]
# bajamos el umbral para las clases más complicadas de predecir.
# Crear una figura y ejes
fig, ax = plt.subplots(figsize=(10, 6))
# Graficar las señales EEG
for i, signal in enumerate(eeg_signals.columns):
ax.plot(eeg_signals.index, eeg_signals[signal] + i * 100, label=f'EEG_{i + 1}')
# Clasificar probabilidades según umbrales
for i, proba_preds in enumerate(proba_predictions):
for j, threshold in enumerate(thresholds):
predicted_class = 1 if proba_preds[j] >= threshold else 0
ax.plot([i, i + 1], [len(eeg_signals.columns) * 100 - (j + 1) * 50, len(eeg_signals.columns) * 100 - (j + 1) * 50],
color='red' if predicted_class == 1 else 'blue', alpha=0.5)
# Ajustar etiquetas y leyenda
ax.set_xlabel('Tiempo')
ax.set_ylabel('EEG Signals')
ax.legend()
# Mostrar la gráfica
plt.show()
Opción (B). Entrenamiento con sólo un experto. Posteriormente, calculamos la curva ROC por experto y luego combiamos los resultados.
from sklearn.preprocessing import MultiLabelBinarizer
from sklearn.multiclass import OneVsRestClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_curve, auc
# Definimos el modelo de Random Forest
classifier = OneVsRestClassifier(
RandomForestClassifier(n_estimators = 1000,
criterion = "gini",
random_state = 0)
)
# Entrenamiento secuencialmente el modelo con los datos de cada sujeto
classifier.fit(combined_X, combined_y_0)
OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))RandomForestClassifier(n_estimators=1000, random_state=0)
RandomForestClassifier(n_estimators=1000, random_state=0)
results_roc_auc_0 = pp.calc_result_roc_auc(classifier, combined_X, combined_y_0, 0)
results_roc_auc_1 = pp.calc_result_roc_auc(classifier, combined_X, combined_y_1, 1)
print ('results_roc_auc_0:',results_roc_auc_0)
print ('results_roc_auc_1:',results_roc_auc_1)
results_roc_auc_0: {0: 1.0, 1: 1.0, 2: 1.0, 3: 1.0, 4: 1.0}
results_roc_auc_1: {0: 0.9918585319159656, 1: 0.9002818508923046, 2: 0.9468489665496386, 3: 0.9839105850723702, 4: 0.9821365881536516}
n_subjects = len(X_train)
n_classes = len(y_train_0[0][0])
print('Datos para el experto 1')
for i in range(n_subjects):
print('\t Subject ', str(i), ';')
for j in range(n_classes):
print('\t\t\t', str(j), ':', results_roc_auc_0[i][j]['roc_auc_train'])
print ('Datos para el experto 2')
for i in range (n_subjects):
print('\t Subject ', str(i), ';')
for j in range(n_classes):
print('\t\t\t', str(j), ':', results_roc_auc_1[i][j]['roc_auc_train'])
Datos para el experto 1 Subject 0 ; 0 : 0.9830332922318124 1 : 0.7531822787141935 2 : 0.7932620416587688 3 : 0.9552282842070838 4 : 0.9527157249233464 Subject 1 ; 0 : 0.9949378151260504 1 : 0.5541295059084275 2 : 0.6540147730895951 3 : 0.9933705955471372 4 : 0.9794366012965037 Subject 2 ; 0 : 0.9781045302551418 1 : 0.817391304347826 2 : 0.7005660865627223 3 : 0.9796797683310842 4 : 0.917712631820266 Subject 3 ; 0 : 0.9849126344086022 1 : 0.7309223623092236 2 : 0.8950975649088856 3 : 0.9675762773532929 4 : 0.8620260468183317 Subject 4 ; 0 : 0.9720426203031765 1 : 0.7429840893703452 2 : 0.8807836803686273 3 : 0.9842302342997755 4 : 0.9298089606212856 Subject 5 ; 0 : 0.9915003501400561 1 : 0.9431292808219178 2 : 0.9384621053019416 3 : 0.9858329725048711 4 : 0.9942604639399772 Datos para el experto 2 Subject 0 ; 0 : 0.9929136561284206 1 : 0.8129262861830541 2 : 0.804018285545561 3 : 0.89724871245164 4 : 0.9430493851007662 Subject 1 ; 0 : 0.9949378151260504 1 : 0.5541295059084275 2 : 0.6540147730895951 3 : 0.9933705955471372 4 : 0.9794366012965037 Subject 2 ; 0 : 0.9805572183976854 1 : 0.7724170369480806 2 : 0.7174747871322357 3 : 0.966895735281379 4 : 0.9099402013151455 Subject 3 ; 0 : 0.9783248247815176 1 : 0.6843520492774224 2 : 0.8957850253721813 3 : 0.9726312821155534 4 : 0.8672727722221536 Subject 4 ; 0 : 0.9645411120382141 1 : 0.6934445069613208 2 : 0.8066956928421897 3 : 0.9880767678579436 4 : 0.9296899645421393 Subject 5 ; 0 : 1.0 1 : 1.0 2 : 1.0 3 : 1.0 4 : 1.0
Matriz de confusión para los datos de entrenamiento (con etiquetas de un experto).
# Datos para el experto 1
print('Datos para el experto 1')
for i in range(n_subjects):
print('\tSubject', str(i), ';')
cm_train_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_train[i], y_train_0[i])
print(cm_train_0)
# Datos para el experto 2
print('Datos para el experto 2')
for i in range(n_subjects):
print('\tSubject', str(i), ';')
cm_train_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_train[i], y_train_1[i])
print(cm_train_1)
Datos para el experto 1 Subject 0 ; [[[613 0] [ 0 301]] [[844 0] [ 0 70]] [[647 0] [ 0 267]] [[719 0] [ 0 195]] [[833 0] [ 0 81]]] Subject 1 ; [[[620 0] [ 0 144]] [[627 0] [ 0 137]] [[530 0] [ 0 234]] [[605 0] [ 0 159]] [[674 0] [ 0 90]]] Subject 2 ; [[[766 0] [ 0 57]] [[690 0] [ 0 133]] [[533 0] [ 0 290]] [[576 0] [ 0 247]] [[727 0] [ 0 96]]] Subject 3 ; [[[714 0] [ 0 80]] [[730 0] [ 0 64]] [[545 0] [ 0 249]] [[496 0] [ 0 298]] [[691 0] [ 0 103]]] Subject 4 ; [[[578 0] [ 0 206]] [[721 0] [ 0 63]] [[629 0] [ 0 155]] [[522 0] [ 0 262]] [[686 0] [ 0 98]]] Subject 5 ; [[[762 0] [ 0 162]] [[810 0] [ 0 114]] [[555 0] [ 0 369]] [[746 0] [ 0 178]] [[823 0] [ 0 101]]] Datos para el experto 2 Subject 0 ; [[[605 26] [ 8 275]] [[784 27] [ 60 43]] [[589 39] [ 58 228]] [[718 51] [ 1 144]] [[814 3] [ 19 78]]] Subject 1 ; [[[617 20] [ 3 124]] [[574 56] [ 53 81]] [[491 55] [ 39 179]] [[586 20] [ 19 139]] [[632 5] [ 42 85]]] Subject 2 ; [[[762 9] [ 4 48]] [[675 82] [ 15 51]] [[454 37] [ 79 253]] [[554 33] [ 22 214]] [[683 3] [ 44 93]]] Subject 3 ; [[[708 13] [ 6 67]] [[701 21] [ 29 43]] [[495 37] [ 50 212]] [[491 34] [ 5 264]] [[673 3] [ 18 100]]] Subject 4 ; [[[577 23] [ 1 183]] [[711 56] [ 10 7]] [[524 12] [105 143]] [[509 39] [ 13 223]] [[673 12] [ 13 86]]] Subject 5 ; [[[754 26] [ 8 136]] [[780 39] [ 30 75]] [[514 27] [ 41 342]] [[724 10] [ 22 168]] [[815 7] [ 8 94]]]
Combinamos las predicciones y comprobmos el auc y matriz de confusión:
# mos los resultados
# combined_results = {}
# for i in range(n_subjects):
# combined_results[i] = {}
# for j in range(n_classes):
# # Combinamos fpr, tpr, thresholds y roc_auc de results_0 y results_1
# combined_results[i][j] = {
# 'fpr': np.concatenate((results_roc_auc_0[i][j]['fpr'], results_roc_auc_1[i][j]['fpr'])),
# 'tpr': np.concatenate((results_roc_auc_0[i][j]['tpr'], results_roc_auc_1[i][j]['tpr'])),
# 'thresholds': np.concatenate((results_roc_auc_0[i][j]['thresholds'], results_roc_auc_1[i][j]['thresholds'])),
# 'roc_auc_train': np.mean([results_roc_auc_0[i][j]['roc_auc_train'], results_roc_auc_1[i][j]['roc_auc_train']])
# }
# for i in range(n_subjects):
# print('Subject ', str(i), ';')
# for j in range(n_classes):
# print('\t\t\t', combined_results[i][j]['roc_auc_train'])
Subject 0 ; 0.9879734741801165 0.7830542824486237 0.7986401636021649 0.9262384983293619 0.9478825550120563 Subject 1 ; 0.9949378151260504 0.5541295059084275 0.6540147730895951 0.9933705955471372 0.9794366012965037 Subject 2 ; 0.9793308743264135 0.7949041706479534 0.7090204368474791 0.9732877518062315 0.9138264165677057 Subject 3 ; 0.9816187295950599 0.707637205793323 0.8954412951405335 0.9701037797344232 0.8646494095202426 Subject 4 ; 0.9682918661706953 0.718214298165833 0.8437396866054085 0.9861535010788596 0.9297494625817124 Subject 5 ; 0.995750175070028 0.9715646404109589 0.9692310526509709 0.9929164862524356 0.9971302319699886
Validación.
# Evaluar el rendimiento en el primer sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu (classifier, X_val [0], y_val_0 [0])
roc_auc_val_1 = pp.calc_roc_acu (classifier, X_val [0], y_val_1 [0])
print("Exactitud para el primer sujeto y primer Experto de validación:", roc_auc_val_0)
print("Exactitud para el primer sujeto y segundo Experto de validación:", roc_auc_val_1)
print ("")
# Realizar la predicción para el segundo sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu (classifier, X_val [0], y_val_0 [0])
roc_auc_val_1 = pp.calc_roc_acu (classifier, X_val [0], y_val_1 [0])
print("Exactitud para el segundo sujeto y primer Exèrto de validación:", roc_auc_val_0)
print("Exactitud para el segundo sujeto y segundo Exèrto de de validación:", roc_auc_val_1)
Exactitud para el primer sujeto y primer Experto de validación: {0: 0.9929254901960785, 1: 0.8497000500804484, 2: 0.9634459028589648, 3: 0.9972391574568117, 4: 0.9847828397191182}
Exactitud para el primer sujeto y segundo Experto de validación: {0: 0.9929254901960785, 1: 0.8497000500804484, 2: 0.9634459028589648, 3: 0.9972391574568117, 4: 0.9847828397191182}
Exactitud para el segundo sujeto y primer Exèrto de validación: {0: 0.9929254901960785, 1: 0.8497000500804484, 2: 0.9634459028589648, 3: 0.9972391574568117, 4: 0.9847828397191182}
Exactitud para el segundo sujeto y segundo Exèrto de de validación: {0: 0.9929254901960785, 1: 0.8497000500804484, 2: 0.9634459028589648, 3: 0.9972391574568117, 4: 0.9847828397191182}
Matriz de confusión para los datos de validación.
for i in range(len(X_val)):
#roc_auc_val_0 = pp.calc_roc_acu(clasificador, X_val[i], y_val_0[i])
#roc_auc_val_1 = pp.calc_roc_acu(clasificador, X_val[i], y_val_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Validación - Sujeto {i + 1}:")
cm_val_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_0[i])
cm_val_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_1[i])
print(f"Experto 1:\n{cm_val_0}")
print(f"Experto 2:\n{cm_val_1}")
Matriz de Confusión para el Conjunto de Validación - Sujeto 1: Experto 1: [[[578 17] [ 12 363]] [[861 0] [109 0]] [[756 20] [ 75 119]] [[807 20] [ 3 140]] [[811 10] [ 30 119]]] Experto 2: [[[578 17] [ 12 363]] [[861 0] [109 0]] [[756 20] [ 75 119]] [[807 20] [ 3 140]] [[811 10] [ 30 119]]] Matriz de Confusión para el Conjunto de Validación - Sujeto 2: Experto 1: [[[592 47] [ 3 124]] [[547 0] [218 1]] [[572 21] [110 63]] [[638 16] [ 10 102]] [[610 21] [ 54 81]]] Experto 2: [[[592 47] [ 3 124]] [[547 0] [218 1]] [[572 21] [110 63]] [[638 16] [ 10 102]] [[610 21] [ 54 81]]]
Conjunto de prueba.
# Evaluar el rendimiento en el primer sujeto de test
roc_auc_test_0 = pp.calc_roc_acu (classifier, X_test [0], y_test_0 [0])
roc_auc_test_1 = pp.calc_roc_acu (classifier, X_test [0], y_test_1 [0])
print("Exactitud para el primer sujeto y primer Experto de Test:", roc_auc_test_0)
print("Exactitud para el primer sujeto y segundo Experto de Test:", roc_auc_test_1)
print ("")
# Realizar la predicción para el segundo sujeto de test
roc_auc_test_0 = pp.calc_roc_acu (classifier, X_test [1], y_test_0 [1])
roc_auc_test_1 = pp.calc_roc_acu (classifier, X_test [1], y_test_1 [1])
print("Exactitud para el segundo sujeto y primer Experto de test:", roc_auc_test_0)
print("Exactitud para el segundo sujeto y segundo Experto de test:", roc_auc_test_1)
Exactitud para el primer sujeto y primer Experto de Test: {0: 0.9805114673836708, 1: 0.823260215320654, 2: 0.8875254938383905, 3: 0.9838530967942732, 4: 0.9535587446546351}
Exactitud para el primer sujeto y segundo Experto de Test: {0: 0.9805114673836708, 1: 0.823260215320654, 2: 0.8875254938383905, 3: 0.9838530967942732, 4: 0.9535587446546351}
Exactitud para el segundo sujeto y primer Experto de test: {0: 0.9889457459926019, 1: 0.8662798194713088, 2: 0.9110275689223057, 3: 0.9678937564873665, 4: 0.9785501533070522}
Exactitud para el segundo sujeto y segundo Experto de test: {0: 0.9941968810278876, 1: 0.8822805578342905, 2: 0.907538833057398, 3: 0.9315350779921637, 4: 0.9721242997875218}
Matriz de confusión par los datos del conjunto de prueba.
for i in range(len(X_test)):
# roc_auc_test_0 = pp.calc_roc_acu(clasificador, X_test[i], y_test_0[i])
# roc_auc_test_1 = pp.calc_roc_acu(clasificador, X_test[i], y_test_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Prueba - Sujeto {i + 1}:")
cm_test_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_test[i], y_test_0[i])
cm_test_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_test[i], y_test_1[i])
print(f"Experto 1:\n{cm_test_0}")
print(f"Experto 2:\n{cm_test_1}")
Matriz de Confusión para el Conjunto de Prueba - Sujeto 1: Experto 1: [[[798 19] [ 30 92]] [[771 1] [166 1]] [[513 64] [150 212]] [[631 83] [ 7 218]] [[862 14] [ 29 34]]] Experto 2: [[[798 19] [ 30 92]] [[771 1] [166 1]] [[513 64] [150 212]] [[631 83] [ 7 218]] [[862 14] [ 29 34]]] Matriz de Confusión para el Conjunto de Prueba - Sujeto 2: Experto 1: [[[796 15] [ 14 86]] [[770 0] [138 3]] [[540 48] [114 209]] [[677 37] [ 35 162]] [[759 2] [ 56 94]]] Experto 2: [[[807 31] [ 3 70]] [[805 0] [103 3]] [[536 53] [118 204]] [[643 25] [ 69 174]] [[744 0] [ 71 96]]]
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [0])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_0 [0].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 9 - Expert 1')
plt.legend(loc='lower right')
plt.show()
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [1])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_1 [1].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 2 - Expert 2')
plt.legend(loc='lower right')
plt.show()
from sklearn.preprocessing import MultiLabelBinarizer
from sklearn.multiclass import OneVsRestClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_curve, auc
n_classes = 5
classifier = OneVsRestClassifier (
RandomForestClassifier(n_estimators = 1000,
criterion = "gini",
random_state = 0)
)
classifier.fit(X_train[0], y_train_0 [0])
rf = []
rf.append(classifier)
y_score_train = classifier.predict_proba(X_train[0])
fpr = dict()
tpr = dict()
thresholds = dict()
roc_auc_train = dict()
for i in range(n_classes):
fpr[i], tpr[i], thresholds[i] = roc_curve(y_train_0[0][:, i], y_score_train [:, i])
roc_auc_train[i] = auc (fpr[i], tpr[i])
roc_auc_train # Con 1 sujeto y los datos de prueba. Con la matriz de características normalizada.
{0: 1.0, 1: 1.0, 2: 1.0, 3: 1.0, 4: 1.0}
# Evaluamos el rendimiento en el primer sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu(classifier, X_val[0], y_val_0[0])
roc_auc_val_1 = pp.calc_roc_acu(classifier, X_val[0], y_val_1[0])
print("Exactitud para el primer sujeto y primer Experto de validación:", roc_auc_val_0)
print("Exactitud para el primer sujeto y segundo Experto de validación:", roc_auc_val_1)
print ("")
# Realizamos la predicción para el segundo sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu(classifier, X_val[0], y_val_0[0])
roc_auc_val_1 = pp.calc_roc_acu(classifier, X_val[0], y_val_1[0])
print("Exactitud para el segundo sujeto y primer Exèrto de validación:", roc_auc_val_0)
print("Exactitud para el segundo sujeto y segundo Exèrto de de validación:", roc_auc_val_1)
Exactitud para el primer sujeto y primer Experto de validación: {0: 0.9797803921568627, 1: 0.5732719581455317, 2: 0.9415287224997343, 3: 0.9963386069794775, 4: 0.9721284405169666}
Exactitud para el primer sujeto y segundo Experto de validación: {0: 0.9797803921568627, 1: 0.5732719581455317, 2: 0.9415287224997343, 3: 0.9963386069794775, 4: 0.9721284405169666}
Exactitud para el segundo sujeto y primer Exèrto de validación: {0: 0.9797803921568627, 1: 0.5732719581455317, 2: 0.9415287224997343, 3: 0.9963386069794775, 4: 0.9721284405169666}
Exactitud para el segundo sujeto y segundo Exèrto de de validación: {0: 0.9797803921568627, 1: 0.5732719581455317, 2: 0.9415287224997343, 3: 0.9963386069794775, 4: 0.9721284405169666}
Matriz de confusión para los datos de validación.
for i in range(len(X_val)):
# roc_auc_val_0 = pp.calc_roc_acu(clasificador, X_val[i], y_val_0[i])
# roc_auc_val_1 = pp.calc_roc_acu(clasificador, X_val[i], y_val_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Validación - Sujeto {i + 1}:")
cm_val_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_0[i])
cm_val_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_1[i])
print(f"Experto 1:\n{cm_val_0}")
print(f"Experto 2:\n{cm_val_1}")
Matriz de Confusión para el Conjunto de Validación - Sujeto 1: Experto 1: [[[560 35] [ 9 366]] [[861 0] [109 0]] [[747 29] [ 70 124]] [[811 16] [ 5 138]] [[815 6] [ 61 88]]] Experto 2: [[[560 35] [ 9 366]] [[861 0] [109 0]] [[747 29] [ 70 124]] [[811 16] [ 5 138]] [[815 6] [ 61 88]]] Matriz de Confusión para el Conjunto de Validación - Sujeto 2: Experto 1: [[[551 88] [ 1 126]] [[547 0] [219 0]] [[562 31] [101 72]] [[651 3] [ 18 94]] [[620 11] [ 64 71]]] Experto 2: [[[551 88] [ 1 126]] [[547 0] [219 0]] [[562 31] [101 72]] [[651 3] [ 18 94]] [[620 11] [ 64 71]]]
# Evaluamos el rendimiento en el primer sujeto de test
roc_auc_test_0 = pp.calc_roc_acu(classifier, X_test[0], y_test_0[0])
roc_auc_test_1 = pp.calc_roc_acu(classifier, X_test[0], y_test_1[0])
print("Exactitud para el primer sujeto y primer Experto de Test:", roc_auc_test_0)
print("Exactitud para el primer sujeto y segundo Experto de Test:", roc_auc_test_1)
print ("")
# Realizamos la predicción para el segundo sujeto de test
roc_auc_test_0 = pp.calc_roc_acu(classifier, X_test[1], y_test_0[1])
roc_auc_test_1 = pp.calc_roc_acu(classifier, X_test[1], y_test_1[1])
print("Exactitud para el segundo sujeto y primer Experto de test:", roc_auc_test_0)
print("Exactitud para el segundo sujeto y segundo Experto de test:", roc_auc_test_1)
Exactitud para el primer sujeto y primer Experto de Test: {0: 0.9649808375303489, 1: 0.7333855604852471, 2: 0.8617300382048507, 3: 0.9614005602240896, 4: 0.9482224396607958}
Exactitud para el primer sujeto y segundo Experto de Test: {0: 0.9649808375303489, 1: 0.7333855604852471, 2: 0.8617300382048507, 3: 0.9614005602240896, 4: 0.9482224396607958}
Exactitud para el segundo sujeto y primer Experto de test: {0: 0.9854315659679409, 1: 0.7908722483190569, 2: 0.8416498178218655, 3: 0.9519046197159066, 4: 0.9757117827420061}
Exactitud para el segundo sujeto y segundo Experto de test: {0: 0.9838166541341093, 1: 0.7973924762686042, 2: 0.8493604277172595, 3: 0.9183669697641752, 4: 0.9619873156911982}
Matriz de confusión para los datos de prueba.
for i in range(len(X_test)):
# roc_auc_test_0 = pp.calc_roc_acu(clasificador, X_test[i], y_test_0[i])
# roc_auc_test_1 = pp.calc_roc_acu(clasificador, X_test[i], y_test_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Prueba - Sujeto {i + 1}:")
cm_test_0 = pp.calculate_multilabel_confusion_matrix(classifier, X_test[i], y_test_0[i])
cm_test_1 = pp.calculate_multilabel_confusion_matrix(classifier, X_test[i], y_test_1[i])
print(f"Experto 1:\n{cm_test_0}")
print(f"Experto 2:\n{cm_test_1}")
Matriz de Confusión para el Conjunto de Prueba - Sujeto 1: Experto 1: [[[758 59] [ 12 110]] [[772 0] [167 0]] [[553 24] [222 140]] [[659 55] [ 18 207]] [[874 2] [ 45 18]]] Experto 2: [[[758 59] [ 12 110]] [[772 0] [167 0]] [[553 24] [222 140]] [[659 55] [ 18 207]] [[874 2] [ 45 18]]] Matriz de Confusión para el Conjunto de Prueba - Sujeto 2: Experto 1: [[[739 72] [ 3 97]] [[770 0] [141 0]] [[558 30] [214 109]] [[709 5] [ 98 99]] [[760 1] [ 86 64]]] Experto 2: [[[742 96] [ 0 73]] [[805 0] [106 0]] [[557 32] [215 107]] [[665 3] [142 101]] [[744 0] [102 65]]]
yasa (yasa.SleepStaging(raw, eeg_name ='C4-A1' , eog_name='LOC-A2', emg_name='X1') (features_2_norm)).from sklearn.preprocessing import label_binarize
def binarize_y(vec):
return label_binarize(vec, classes = [0, 1, 2, 3, 4])
X_train = [sub.get_yasa_optim_features().values for sub in X_train_o]
X_val = [sub.get_yasa_optim_features().values for sub in X_val_o]
X_test = [sub.get_yasa_optim_features().values for sub in X_test_o]
y_train_0 = [binarize_y(sub.get_labels()[0]) for sub in X_train_o]
y_train_1 = [binarize_y(sub.get_labels()[1]) for sub in X_train_o]
y_val_0 = [binarize_y(sub.get_labels()[0]) for sub in X_val_o]
y_val_1 = [binarize_y(sub.get_labels()[1]) for sub in X_val_o]
y_test_0 = [binarize_y(sub.get_labels()[0]) for sub in X_test_o]
y_test_1 = [binarize_y(sub.get_labels()[1]) for sub in X_test_o]
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import label_binarize
from sklearn.multiclass import OneVsRestClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_curve, auc
# Combinamos las características y etiquetas de los seis sujetos en una sola estructura
combined_X = []
combined_y_0 = []
combined_y_1 = []
for i in range(6):
combined_X.extend(X_train[i])
combined_y_0.extend(y_train_0[i])
combined_y_1.extend(y_train_1[i])
combined_X = np.array(combined_X)
combined_y_0 = np.array(combined_y_0)
combined_y_1 = np.array(combined_y_1)
combined_x_all = np.concatenate([combined_X, combined_X])
combined_y_all = np.concatenate([combined_y_0, combined_y_1])
classifier = OneVsRestClassifier(
RandomForestClassifier(n_estimators = 1000,
criterion = "gini",
random_state = 0)
)
classifier.fit(combined_x_all, combined_y_all)
OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))RandomForestClassifier(n_estimators=1000, random_state=0)
RandomForestClassifier(n_estimators=1000, random_state=0)
from sklearn.inspection import permutation_importance
for i in range (len(X_test)):
for j in range (len(y_test_0)):
# Calcular la importancia de la permutación
result = permutation_importance(classifier, X_test [i], y_test_0 [j], n_repeats = 2, random_state = 42)
# Obtener las importancias
importances = result.importances_mean
# Obtener índices ordenados por importancia
indices = importances.argsort()[::-1]
print (f'Las 30 features más importantes para el subject {i} del Experto 1:')
# Imprimir las características ordenadas por importancia
for f in range(X_test[i].shape [1]):
print(f"\t\t{X_test[i].columns[indices[f]]}: {importances[indices[f]]}")
# Como tarda bastante, la traza que saccamos es bajo el subject 0 y experto 1 . no se hace el for...
# si se desea para los dos subjects hay que esperar bastante y consumir muchisimos recursos ciputacionales...
from sklearn.inspection import permutation_importance
for i in range (len(X_test)):
for j in range (len(y_test_0)):
# Calcular la importancia de la permutación
result = permutation_importance(classifier, X_test [i], y_test_1 [j], n_repeats = 30, random_state = 42)
# Obtener las importancias
importances = result.importances_mean
# Obtener índices ordenados por importancia
indices = importances.argsort()[::-1]
print (f' Las 30 features más importantes para el subject {i} del Experto 2:)
# Imprimir las características ordenadas por importancia
for f in range(X_test[i] .shape[1]):
print(f"{X_test[i].columns[indices[f]]}: {importances[indices[f]]}")
pp.calc_roc_acu (classifier, combined_x_all, combined_y_all)
{0: 0.9992726157134799,
1: 0.9881299416149647,
2: 0.9924173893957494,
3: 0.9981082199322366,
4: 0.9985678162262597}
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [0])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_0 [0].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 9 - Expert 1')
plt.legend(loc='lower right')
plt.show()
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [1])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_1 [1].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 2 - Expert 2')
plt.legend(loc='lower right')
plt.show()
Matriz de confusión para los datos de entrenamiento.
cm_entrenamiento = pp.calculate_multilabel_confusion_matrix (classifier, combined_x_all, combined_y_all)
print("Matriz de Confusión para el Conjunto de Entrenamiento:")
print(cm_entrenamiento)
Matriz de Confusión para el Conjunto de Entrenamiento: [[[8125 68] [ 79 1734]] [[8857 71] [ 407 671]] [[6438 275] [ 304 2989]] [[7326 107] [ 162 2411]] [[8653 104] [ 73 1176]]]
Validación.
# Evaluar el rendimiento en el primer sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu (classifier, X_val [0], y_val_0 [0])
roc_auc_val_1 = pp.calc_roc_acu (classifier, X_val [0], y_val_1 [0])
print("Exactitud para el primer sujeto y primer Experto de validación:", roc_auc_val_0)
print("Exactitud para el primer sujeto y segundo Experto de validación:", roc_auc_val_1)
print ("")
# Realizar la predicción para el segundo sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu (classifier, X_val [1], y_val_0 [1])
roc_auc_val_1 = pp.calc_roc_acu (classifier, X_val [1], y_val_1 [1])
print("Exactitud para el segundo sujeto y primer Experto de validación:", roc_auc_val_0)
print("Exactitud para el segundo sujeto y segundo Experto de validación:", roc_auc_val_1)
Exactitud para el primer sujeto y primer Experto de validación: {0: 0.9925557422969187, 1: 0.8440260418331577, 2: 0.9458165054734828, 3: 0.9985540457124495, 4: 0.99664020796377}
Exactitud para el primer sujeto y segundo Experto de validación: {0: 0.9925557422969187, 1: 0.8440260418331577, 2: 0.9458165054734828, 3: 0.9985540457124495, 4: 0.99664020796377}
Exactitud para el segundo sujeto y primer Experto de validación: {0: 0.9660024891254791, 1: 0.89774444249664, 2: 0.8902416438409576, 3: 0.9976381607688947, 4: 0.9803603920878089}
Exactitud para el segundo sujeto y segundo Experto de validación: {0: 0.9660024891254791, 1: 0.89774444249664, 2: 0.8902416438409576, 3: 0.9976381607688947, 4: 0.9803603920878089}
Matriz de confusión para los datos de validación:
for i in range(len(X_val)):
#roc_auc_val_0 = pp.calc_roc_acu(clasificador, X_val[i], y_val_0[i])
#roc_auc_val_1 = pp.calc_roc_acu(clasificador, X_val[i], y_val_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Validación - Sujeto {i + 1}:")
cm_val_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_0[i])
cm_val_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_1[i])
print(f"Experto 1:\n{cm_val_0}")
print(f"Experto 2:\n{cm_val_1}")
Matriz de Confusión para el Conjunto de Validación - Sujeto 1: Experto 1: [[[577 18] [ 37 338]] [[861 0] [107 2]] [[703 73] [ 39 155]] [[812 15] [ 2 141]] [[812 9] [ 13 136]]] Experto 2: [[[577 18] [ 37 338]] [[861 0] [107 2]] [[703 73] [ 39 155]] [[812 15] [ 2 141]] [[812 9] [ 13 136]]] Matriz de Confusión para el Conjunto de Validación - Sujeto 2: Experto 1: [[[570 69] [ 12 115]] [[546 1] [191 28]] [[562 31] [ 89 84]] [[650 4] [ 9 103]] [[631 0] [ 83 52]]] Experto 2: [[[570 69] [ 12 115]] [[546 1] [191 28]] [[562 31] [ 89 84]] [[650 4] [ 9 103]] [[631 0] [ 83 52]]]
Test.
# Evaluar el rendimiento en el primer sujeto de test
roc_auc_test_0 = pp.calc_roc_acu (classifier, X_test [0], y_test_0 [0])
roc_auc_test_1 = pp.calc_roc_acu (classifier, X_test [0], y_test_1 [0])
print("Exactitud para el primer sujeto y primer Experto de Test:", roc_auc_test_0)
print("Exactitud para el primer sujeto y segundo Experto de Test:", roc_auc_test_1)
print ("")
# Realizar la predicción para el segundo sujeto de test
roc_auc_test_0 = pp.calc_roc_acu (classifier, X_test [1], y_test_0 [1])
roc_auc_test_1 = pp.calc_roc_acu (classifier, X_test [1], y_test_1 [1])
print("Exactitud para el segundo sujeto y primer Experto de test:", roc_auc_test_0)
print("Exactitud para el segundo sujeto y segundo Experto de test:", roc_auc_test_1)
Exactitud para el primer sujeto y primer Experto de Test: {0: 0.9730421173024059, 1: 0.8654982780552883, 2: 0.9304461062650211, 3: 0.9899844382197324, 4: 0.9862832499818802}
Exactitud para el primer sujeto y segundo Experto de Test: {0: 0.9730421173024059, 1: 0.8654982780552883, 2: 0.9304461062650211, 3: 0.9899844382197324, 4: 0.9862832499818802}
Exactitud para el segundo sujeto y primer Experto de test: {0: 0.9858199753390876, 1: 0.8600027631942525, 2: 0.93078283945157, 3: 0.9895562285828037, 4: 0.9890757774857644}
Exactitud para el segundo sujeto y segundo Experto de test: {0: 0.9946055513780365, 1: 0.9080042189148012, 2: 0.9227715150428667, 3: 0.9681501195140583, 4: 0.9899756937737428}
Matriz de confusión para los datos de prueba.
for i in range(len(X_test)):
#roc_auc_test_0 = pp.calc_roc_acu(clasificador, X_test[i], y_test_0[i])
#roc_auc_test_1 = pp.calc_roc_acu(clasificador, X_test[i], y_test_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Prueba - Sujeto {i + 1}:")
cm_test_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_test[i], y_test_0[i])
cm_test_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_test[i], y_test_1[i])
print(f"Experto 1:\n{cm_test_0}")
print(f"Experto 2:\n{cm_test_1}")
Matriz de Confusión para el Conjunto de Prueba - Sujeto 1: Experto 1: [[[795 22] [ 38 84]] [[770 2] [159 8]] [[481 96] [ 28 334]] [[704 10] [ 37 188]] [[861 15] [ 12 51]]] Experto 2: [[[795 22] [ 38 84]] [[770 2] [159 8]] [[481 96] [ 28 334]] [[704 10] [ 37 188]] [[861 15] [ 12 51]]] Matriz de Confusión para el Conjunto de Prueba - Sujeto 2: Experto 1: [[[791 20] [ 11 89]] [[770 0] [138 3]] [[526 62] [ 62 261]] [[708 6] [ 28 169]] [[751 10] [ 41 109]]] Experto 2: [[[800 38] [ 2 71]] [[804 1] [104 2]] [[522 67] [ 66 256]] [[666 2] [ 70 173]] [[740 4] [ 52 115]]]
from sklearn.preprocessing import MultiLabelBinarizer
from sklearn.multiclass import OneVsRestClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_curve, auc
# Definimos el modelo de Random Forest
classifier = OneVsRestClassifier(
RandomForestClassifier(n_estimators = 1000,
criterion = "gini",
random_state = 0)
)
# Entrenamos secuencialmente el modelo con los datos de cada sujeto
classifier.fit(combined_X, combined_y_0)
OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. OneVsRestClassifier(estimator=RandomForestClassifier(n_estimators=1000,
random_state=0))RandomForestClassifier(n_estimators=1000, random_state=0)
RandomForestClassifier(n_estimators=1000, random_state=0)
results_roc_auc_0 = pp.calc_roc_acu (classifier, combined_X, combined_y_0)
results_roc_auc_1 = pp.calc_roc_acu (classifier, combined_X, combined_y_1)
print ('results_roc_auc_0:',results_roc_auc_0)
print ('results_roc_auc_1:',results_roc_auc_1)
results_roc_auc_0: {0: 1.0, 1: 0.9999999999999999, 2: 1.0, 3: 1.0, 4: 1.0000000000000002}
results_roc_auc_1: {0: 0.9934025503663773, 1: 0.9107088156993448, 2: 0.9524624492955522, 3: 0.9892074644599185, 4: 0.9933002000244928}
# Datos para el experto 1
print('Datos para el experto 1')
for i in range(n_subjects):
print('\tSubject', str(i), ';')
cm_train_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_train[i], y_train_0[i])
print(cm_train_0)
# Datos para el experto 2
print('Datos para el experto 2')
for i in range(n_subjects):
print('\tSubject', str(i), ';')
cm_train_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_train[i], y_train_1[i])
print(cm_train_1)
Datos para el experto 1 Subject 0 ; [[[613 0] [ 0 301]] [[844 0] [ 0 70]] [[647 0] [ 0 267]] [[719 0] [ 0 195]] [[833 0] [ 0 81]]] Subject 1 ; [[[620 0] [ 0 144]] [[627 0] [ 0 137]] [[530 0] [ 0 234]] [[605 0] [ 0 159]] [[674 0] [ 0 90]]] Subject 2 ; [[[766 0] [ 0 57]] [[690 0] [ 0 133]] [[533 0] [ 0 290]] [[576 0] [ 0 247]] [[727 0] [ 0 96]]] Subject 3 ; [[[714 0] [ 0 80]] [[730 0] [ 0 64]] [[545 0] [ 0 249]] [[496 0] [ 0 298]] [[691 0] [ 0 103]]] Subject 4 ; [[[578 0] [ 0 206]] [[721 0] [ 0 63]] [[629 0] [ 0 155]] [[522 0] [ 0 262]] [[686 0] [ 0 98]]] Subject 5 ; [[[762 0] [ 0 162]] [[810 0] [ 0 114]] [[555 0] [ 0 369]] [[746 0] [ 0 178]] [[823 0] [ 0 101]]] Datos para el experto 2 Subject 0 ; [[[605 26] [ 8 275]] [[784 27] [ 60 43]] [[589 39] [ 58 228]] [[718 51] [ 1 144]] [[814 3] [ 19 78]]] Subject 1 ; [[[617 20] [ 3 124]] [[574 56] [ 53 81]] [[491 55] [ 39 179]] [[586 20] [ 19 139]] [[632 5] [ 42 85]]] Subject 2 ; [[[762 9] [ 4 48]] [[675 82] [ 15 51]] [[454 37] [ 79 253]] [[554 33] [ 22 214]] [[683 3] [ 44 93]]] Subject 3 ; [[[708 13] [ 6 67]] [[701 21] [ 29 43]] [[495 37] [ 50 212]] [[491 34] [ 5 264]] [[673 3] [ 18 100]]] Subject 4 ; [[[577 23] [ 1 183]] [[711 56] [ 10 7]] [[524 12] [105 143]] [[509 39] [ 13 223]] [[673 12] [ 13 86]]] Subject 5 ; [[[754 26] [ 8 136]] [[780 39] [ 30 75]] [[514 27] [ 41 342]] [[724 10] [ 22 168]] [[815 7] [ 8 94]]]
# Evaluar el rendimiento en el primer sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu (classifier, X_val [0], y_val_0 [0])
roc_auc_val_1 = pp.calc_roc_acu (classifier, X_val [0], y_val_1 [0])
print("Exactitud para el primer sujeto y primer Experto de validación:", roc_auc_val_0)
print("Exactitud para el primer sujeto y segundo Experto de validación:", roc_auc_val_1)
print ("")
# Realizar la predicción para el segundo sujeto de validación
roc_auc_val_0 = pp.calc_roc_acu (classifier, X_val [1], y_val_0 [1])
roc_auc_val_1 = pp.calc_roc_acu (classifier, X_val [1], y_val_1 [1])
print("Exactitud para el segundo sujeto y primer Experto de validación:", roc_auc_val_0)
print("Exactitud para el segundo sujeto y segundo Experto de validación:", roc_auc_val_1)
Exactitud para el primer sujeto y primer Experto de validación: {0: 0.9924907563025209, 1: 0.8269667231403638, 2: 0.9323387182484855, 3: 0.9981819872992788, 4: 0.9974781123037055}
Exactitud para el primer sujeto y segundo Experto de validación: {0: 0.9924907563025209, 1: 0.8269667231403638, 2: 0.9323387182484855, 3: 0.9981819872992788, 4: 0.9974781123037055}
Exactitud para el segundo sujeto y primer Experto de validación: {0: 0.9646408635540276, 1: 0.894355262828379, 2: 0.8978691672596477, 3: 0.997924858016601, 4: 0.9834653988378236}
Exactitud para el segundo sujeto y segundo Experto de validación: {0: 0.9646408635540276, 1: 0.894355262828379, 2: 0.8978691672596477, 3: 0.997924858016601, 4: 0.9834653988378236}
Matriz de confusión para los datos de validación:
for i in range(len(X_val)):
#roc_auc_val_0 = pp.calc_roc_acu(clasificador, X_val[i], y_val_0[i])
#roc_auc_val_1 = pp.calc_roc_acu(clasificador, X_val[i], y_val_1[i])
print(f"\nMatriz de Confusión para el Conjunto de Validación - Sujeto {i + 1}:")
cm_val_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_0[i])
cm_val_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_val[i], y_val_1[i])
print(f"Experto 1:\n{cm_val_0}")
print(f"Experto 2:\n{cm_val_1}")
Matriz de Confusión para el Conjunto de Validación - Sujeto 1: Experto 1: [[[577 18] [ 35 340]] [[861 0] [108 1]] [[705 71] [ 52 142]] [[787 40] [ 0 143]] [[813 8] [ 12 137]]] Experto 2: [[[577 18] [ 35 340]] [[861 0] [108 1]] [[705 71] [ 52 142]] [[787 40] [ 0 143]] [[813 8] [ 12 137]]] Matriz de Confusión para el Conjunto de Validación - Sujeto 2: Experto 1: [[[561 78] [ 11 116]] [[547 0] [213 6]] [[568 25] [ 89 84]] [[648 6] [ 4 108]] [[631 0] [ 86 49]]] Experto 2: [[[561 78] [ 11 116]] [[547 0] [213 6]] [[568 25] [ 89 84]] [[648 6] [ 4 108]] [[631 0] [ 86 49]]]
Conjunto de prueba.
# Evaluar el rendimiento en el primer sujeto de test
roc_auc_test_0 = pp.calc_roc_acu (classifier, X_test [0], y_test_0 [0])
roc_auc_test_1 = pp.calc_roc_acu (classifier, X_test [0], y_test_1 [0])
print("Exactitud para el primer sujeto y primer Experto de Test:", roc_auc_test_0)
print("Exactitud para el primer sujeto y segundo Experto de Test:", roc_auc_test_1)
print ("")
# Realizar la predicción para el segundo sujeto de test
roc_auc_test_0 = pp.calc_roc_acu (classifier, X_test [1], y_test_0 [1])
roc_auc_test_1 = pp.calc_roc_acu (classifier, X_test [1], y_test_1 [1])
print("Exactitud para el segundo sujeto y primer Experto de test:", roc_auc_test_0)
print("Exactitud para el segundo sujeto y segundo Experto de test:", roc_auc_test_1)
Exactitud para el primer sujeto y primer Experto de Test: {0: 0.9725053674980436, 1: 0.8525759362104804, 2: 0.9190181640606297, 3: 0.9899719887955183, 4: 0.9882673769660071}
Exactitud para el primer sujeto y segundo Experto de Test: {0: 0.9725053674980436, 1: 0.8525759362104804, 2: 0.9190181640606297, 3: 0.9899719887955183, 4: 0.9882673769660071}
Exactitud para el segundo sujeto y primer Experto de test: {0: 0.9853884093711467, 1: 0.8443768996960486, 2: 0.9294902171394874, 3: 0.9870217122381948, 4: 0.9888961892247043}
Exactitud para el segundo sujeto y segundo Experto de test: {0: 0.9938372511197568, 1: 0.8829426930739482, 2: 0.9217011673644139, 3: 0.9582440058155294, 4: 0.9906598094134311}
Matriz de confusión para los datos de prueba.
for i in range(len(X_test)):
#roc_auc_test_0 = pp.calc_roc_acu(clasificador, X_test[i], y_test_0[i])
#roc_auc_test_1 = pp.calc_roc_acu(clasificador, X_test[i], y_test_1[i])
print(f"\nMatriz de Confusión para el Conjunto de test - Sujeto {i + 1}:")
cm_test_0 = pp.calculate_multilabel_confusion_matrix (classifier, X_test[i], y_test_0[i])
cm_test_1 = pp.calculate_multilabel_confusion_matrix (classifier, X_test[i], y_test_1[i])
print(f"Experto 1:\n{cm_test_0}")
print(f"Experto 2:\n{cm_test_1}")
Matriz de Confusión para el Conjunto de test - Sujeto 1: Experto 1: [[[792 25] [ 32 90]] [[771 1] [161 6]] [[483 94] [ 29 333]] [[704 10] [ 31 194]] [[861 15] [ 15 48]]] Experto 2: [[[792 25] [ 32 90]] [[771 1] [161 6]] [[483 94] [ 29 333]] [[704 10] [ 31 194]] [[861 15] [ 15 48]]] Matriz de Confusión para el Conjunto de test - Sujeto 2: Experto 1: [[[787 24] [ 11 89]] [[770 0] [141 0]] [[528 60] [ 64 259]] [[701 13] [ 33 164]] [[755 6] [ 45 105]]] Experto 2: [[[796 42] [ 2 71]] [[805 0] [106 0]] [[521 68] [ 71 251]] [[661 7] [ 73 170]] [[742 2] [ 58 109]]]
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [0])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_0 [0].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 9 - Expert 1')
plt.legend(loc='lower right')
plt.show()
# Calcula las probabilidades de las predicciones en el conjunto de prueba
y_prob = classifier.predict_proba(X_test [1])
# Calcula la curva ROC
fpr, tpr, _ = roc_curve(y_test_1 [1].ravel(), y_prob.ravel())
roc_auc = auc(fpr, tpr)
# Grafica la curva ROC
plt.figure(figsize=(8, 8))
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='Random Guess')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve Subkject 2 - Expert 2')
plt.legend(loc='lower right')
plt.show()
En la evaluación de nuestro modelo para la identificación de etapas del sueño, se han empleado dos estrategias de entrenamiento con diferentes conjuntos de características, así como la división de datos en conjuntos de entrenamiento, validación y prueba utilizando la estrategia hold-out. Se destacan dos estrategias de entrenamiento: una que utiliza las etiquetas de dos expertos y otra que utiliza solo las etiquetas de un experto. Además, se ha empleado la librería YASA para extraer características relevantes.
En base a los resultados anteirores, utilizando las etiquetas de dos expertos, se muestra un rendimiento impresionante en la métrica del área bajo la curva (AUC) en el conjunto de entrenamiento y validación. Sin embargo, al analizar las matrices de confusión y métricas específicas como sensibilidad y especificidad, se revelan aspectos importantes.
Se destaca que el modelo presenta una alta especificidad en la fase N1, indicando que es efectivo para predecir correctamente esta etapa del sueño. Sin embargo, la sensibilidad en la fase N2 disminuye, lo que sugiere que el modelo tiene dificultades para predecir esta etapa, pero mantiene una alta especificidad en otras fases. Se resalta que, en el contexto del sueño, una baja sensibilidad en la fase N2 puede no ser crítica para la salud del sujeto, ya que esta etapa corresponde a un sueño ligero.
El análisis de la estrategia de utilizar solo las etiquetas de un experto muestra resultados más desafiantes, especialmente en la predicción de la fase N2 y REM. Las métricas de sensibilidad y AUC disminuyen en comparación con la estrategia anterior. Se observa que la sensibilidad en la fase N2 disminuye aún más, lo que puede deberse a las variaciones en las predicciones de los expertos utilizados para el entrenamiento.
La evaluación con la librería YASA muestra resultados comparables, aunque ligeramente mejores en algunas métricas. Sin embargo, se destaca que la sensibilidad puede variar según el sujeto, y el modelo puede tener dificultades, especialmente en la fase N2.
Finalmente, se identifica al Subject 8 como el más desafiante de predecir, con bajos niveles de sensibilidad en varias etapas del sueño. Se resalta la importancia de mantener una alta especificidad para garantizar que, cuando el modelo predice una etapa, sea altamente probable que sea correcta.
En resumen, el modelo demuestra ser eficaz en la identificación de ciertas etapas del sueño, con fortalezas en la fase N1 y limitaciones en la sensibilidad de la fase N2. Se reconoce la importancia de mantener una alta especificidad, especialmente en el contexto de la clasificación del sueño, donde algunas etapas pueden no tener un impacto crítico en la salud.
La visión clínica de la solución planteada se centra en evaluar la utilidad y aplicabilidad del modelo de clasificación en el contexto de la medicina del sueño. Algunos aspectos clave incluyen:
Especificidad y Sensibilidad:
Contexto Clínico:
Resultados Detallados y Optimización:
Dificultades en Predicción de Etapas Específicas:
Interpretabilidad del Modelo:
Validación Externa:
En resumen, la visión clínica de la solución destaca la importancia de evaluar la eficiencia y robustez del modelo desde una perspectiva práctica y clínica. La combinación de alta especificidad, sensibilidad aceptable y una fácil interpretabilidad del modelo respalda su idoneidad para aplicaciones médicas relacionadas con la identificación de etapas del sueño.
La selección de variables basada en las marcas de los expertos sobre las fases de sueño es una estrategia de interpretabilidad que implica utilizar únicamente aquellas características que los expertos consideran relevantes para la identificación de las diferentes etapas del sueño. Este enfoque tiene varias ventajas en términos de explicabilidad del modelo:
Relevancia Clínica Directa:
Simplicidad y Claridad:
Facilita la Comparación con Decisiones Expertas:
Reducción del Riesgo de Overfitting:
Explicabilidad Directa:
Aunque este enfoque simplifica el modelo y mejora su explicabilidad, es esencial considerar que también puede haber información valiosa en otras características que no fueron marcadas por expertos. Por lo tanto, es recomendable equilibrar la simplificación del modelo con la necesidad de capturar la complejidad inherente de los datos. Además, la iteración continua con expertos en medicina del sueño puede ser valiosa para ajustar y mejorar el modelo en base a nuevos conocimientos y descubrimientos clínicos.
La solución planteada demuestra ser prometedora desde una perspectiva clínica. Se destaca la alta especificidad del modelo, lo que indica su capacidad para predecir correctamente las etapas de sueño que no pertenecen a la clase de interés, en este caso, la fase "No REM". Además, se observa una sensibilidad aceptable en la mayoría de las etapas, lo que sugiere que el modelo tiene una baja tendencia a clasificar erróneamente las etapas de sueño.
La robustez del modelo se verifica al combinar todas las variables y etiquetas de los dos expertos en una única estructura, lo que resulta en un rendimiento mejorado en términos de especificidad y sensibilidad en la mayoría de las etapas. Sin embargo, se reconoce que hay subjects que son más difíciles de predecir, y los resultados de la curva ROC son menores en esos casos. Se atribuye esta variabilidad a posibles interferencias, ruido en las señales o características particulares de los hábitos de sueño de esos individuos.
La visión clínica de la solución planteada enfatiza la importancia de considerar el contexto clínico al interpretar la eficiencia y robustez del modelo. La detección precisa de la fase REM puede ser crítica en ciertos contextos clínicos, y se sugiere evaluar la eficiencia del modelo en función de la relevancia clínica de cada clase. La interpretabilidad y explicabilidad del modelo se destacan como aspectos positivos, ya que se emplea un algoritmo de selección automática de características que es fácilmente explicable.
La conclusión general es que el modelo cumple con los requisitos clínicos y demuestra ser robusto en diversas situaciones. Aunque se señala que la aplicación puede ser mejorable, con sugerencias como la eliminación de artefactos y la comparación con otros modelos como el Gradient Boosting de Yasa, se destaca que el modelo es capaz de predecir correctamente las etapas de sueño con un área bajo la curva (AUC) significativo. Además, su fácil explicabilidad es un punto a favor, y se sugiere la posibilidad de realizar validaciones externas con expertos en medicina del sueño para garantizar su aplicabilidad y precisión en entornos clínicos reales. En general, la solución presenta una base sólida para la identificación de etapas del sueño, con potencial para mejoras adicionales y aplicaciones clínicas significativas.